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ABSTRACT 

A  discussion  of  the  analytical  methods  of  pre¬ 
dicting  uniaxial  stress  wave  propagation  and  spall  as  a  re¬ 
sult  of  uniformly  distributed  intense  impulsive  loading  is 
given.  For  the  materials  for  which  data  currently  exists, 
no  strain-rate  effects  are  observed  on  yield  (i.e.  two 
aluminum  alloys)  and  the  strain-rate  effect  on  the  fracture 
threshold  is  found  to  be  small  (copper) .  For  such  materials 
a  detailed  analytical  treatment  is  possible,  once  the  con¬ 
stants  associated  with  the  constitutive  equation  and  frac¬ 
ture  criterion  are  evaluated.  The  plate  impact  experiment 
provides  a  convenient  means  of  evaluating  these  constants. 
Considerations  entering  into  the  interpretation  of  the  plate 
impact  experiment  and  a  detailed  analysis  of  the  effect  of 
air  between  the  impacting  plates  are  given. 


This  report  has  been  reviewed  and  is  approved. 
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SECTION  I 


INTRODUCTION 


The  prediction  of  the  response  of  shell- like  struc¬ 
tures  to  intense  distributed  impulsive  loads  involves  sever¬ 
al  distinct  problems.  In  the  initial  phase,  stress  waves 
are  propagated  through  the  thickness  of  the  structure.  Re¬ 
flection  of  these  stress  waves  at  an  inner  free  surface, 
and  interaction  of  stress  waves,  may  lead  to  transient  ten¬ 
sile  stresses  of  sufficient  magnitude  to  induce  fracture  and 
lead  to  spallation.  In  the  case  of  multilayer  structures, 
transient  tensile  stresses  may  lead  to  unbonding  and  sepa¬ 
ration  of  layers. 

The  subsequent  response  of  the  structure  may  in¬ 
clude  elastic  vibration,  or,  in  more  severe  cases,  large 
permanent  plastic  deformation,  possibly  accompanied  by 
failure  of  supporting  substructure  and/or  tearing  of  the 
skin. 

In  predicting  the  total  response  of  the  structure, 
it  is  first  necessary  to  carry  out  an  analysis  of  stress 
wave  propagation  in  order  to  determine  if  spall  occurs,  and 
to  determine  the  kinetic  energy  distribution  among  the  spall 
and  remaining  separated  material  layers.  This  information 
can  then  be  utilized  in  an  analysis  of  the  subsequent  deforma¬ 
tion  and  rupture  of  the  structure. 

The  present  report  is  concerned  only  with  the  stress 
wave  propagation  and  spallation  problem.  The  deformation  prob¬ 
lem  is  considered  in  another  report  generated  under  the  same 
1* 

contract . 


Superscripts  refer  to  references  listed  at  the  end  of  the 
report. 

Manuscript  released  by  the  authors  April  1962  for  publication 
as  an  ASD  Technical  Documentary  Report. 
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The  impulsive  loads  considered  in  this  study  are, 
in  general, distributed  over  the  surface  of  the  structure  in 
such  a  way  that  the  variation  in  intensity  over  an  area  of 
the  surface  with  dimensions  corresponding  to  several  times 
the  thickness  of  the  structure  is  negligible.  Thus  the 
problem  may  be  treated  as  locally  one - dimen sional,  with 
slab  symmetry  if  the  structure  is  locally  nearly  flat,  or 
with  cylindrical  or  spherical  symmetry  if  the  structure  is 
locally  curved.  The  solution  is  then  built  up  of  a  number 
of  locally  one -dimensional  calculations.  Methods  of  solu¬ 
tion,  in  cases  where  local  one- dimensional  symmetry  cannot 
be  assumed,  generally  become  so  unwieldy  that  the  above 
method  leads  to  great  simplification  when  it  can  be  used. 

Loads  to  be  considered  induce  very  intense 
stresses  of  very  short  duration.  The  "equation  of  state" 
of  the  structural  materials  in  this  range  is  strongly  non¬ 
linear,  necessitating  numerical  methods  of  solution.  On 
the  other  hand,  stresses  are  found  to  attenuate  rapidly  to 
levels  where  elastic  effects  become  appreciable,  and  elastic- 
plastic  transitions  must  be  considered. 

The  wave  propagation  problem  has  been  handled  by 
two  methods,  each  of  which  has  special  merits  making  it  use¬ 
ful  in  specific  instances.  The  two  methods  are;  a  straight 
numerical  solution  of  the  nonlinear  differential  equations 
of  motion  in  finite  difference  form,  and  a  characteristics 
method.  While  both  methods  are  well  known,  a  section  is 
devoted  to  each  in  the  following,  outlining  details  which 
are  peculiar  to  the  methods  used  here. 

The  equation  of  state  of  structural  materials  is 
discussed  in  some  detail,  particularly  in  regard  to  exten¬ 
sion  to  elastic-plastic  behavior.  Copper  and  aluminum  are 
considered  as  quantitative  examples.  Consequences  of  the 
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equations  of  state,  and  of  the  use  of  the  two  methods  of 
solution  are  investigated. 

Using  the  above  analytical  methods,  spall  frac¬ 
ture  criteria  are  examined,  with  a  quantitative  example  for 

22 

copper  based  on  data  of  Lundergan,  and  the  effect  of  air 
between  impacting  flat  plates  is  analyzed. 
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SECTION  II 


FINITE  DIFFERENCE  METHOD 


It  is  convenient  to  write  the  equations  expressing 
conservation  of  mass,  momentum*  and  energy  in  a  one-dimen¬ 
sional  medium  in  Lagrangian  form.  The  coordinate  system  is 
thereby  fixed  in  the  medium  and  distorts  with  it.  Problems 
associated  with  the  transport  of  mass  relative  to  the  coordi¬ 
nate  system  are  thereby  eliminated.  Consequently*  the 
(Eulerian)  position  R  of  a  particle  is  taken  to  be  a  func¬ 
tion  of  time  and  Lagrangian  coordinate  r*  here  taken  to  be 
the  initial  position  of  the  particle  at  time  zero.  It  is 
generally  not  convenient  to  eliminate  R  altogether  in  a 
finite  difference  scheme  and  find  the  solution  solely  in 
terms  of  r  and  t.  Consequently*  the  Eulerian  coordinate 
R(r,t)  will  be  retained  throughout  the  formulation. 

The  differential  equations  of  motion  are  well 
known.  They  are  derived  by  considering  an  element  of  the 
medium*  and  applying  conservation  principles.  Generalizing 
to  slab,  cylindrical  and  spherical  symmetry,  the  element 
becomes  a  flat  slice,  taken  to  be  of  lateral  area  tt,  a 
cylindrical  shell  of  unit  lateral  width*  with  surface  area 
2ttR*  or  a  spherical  shell  of  surface  area  4ttR2*  respectively. 
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The  surface  area  of  the  element  at  position  R 
(2-1 

thus  becomes  tt(2R)  where  a  =  1  for  slab,  2  for  cylindri¬ 
cal  and  3  for  spherical  symmetry. 

Conservation  of  mass  is  expressed  by  the  state¬ 
ment  that  the  time  rate  of  change  of  mass  in  a  Lagrangian 
element  is  zero,  i.e.. 


I if  V’W'lfdr-O  (2.1a) 

ri 

where  specific  volume  is  denoted  by  V. 

Conservation  of  momentum  is  expressed  by  the 
statement  that  the  time  rate  of  change  of  (radial)  momentum 
of  the  element  equals  the  net  force  on  the  element,  i.e., 

Vz 

a  r  i  a_1  ad  a_i 

It  J  v  u  lr(2R)  H  dr  =  (Pl  +  Ql)  ir<2Ri) 

ri  (2.1b) 

a-i 

-  (P2  +  Qs>)  tt(2R2) 

where  U  denotes  velocity,  P  denotes  pressure,  and  Q  is  an 
additional  term  denoting  viscous  forces. 

Conservation  of  energy  is  expressed  in  the  state¬ 
ment  that  the  time  rate  of  change  of  kinetic  and  internal 
energy  of  the  segment  must  equal  the  net  power  flowing  into 
the  segment,  i.e.. 

It/  V<7  U*  +  ®  "  W)r(2R)“  1  H  dr 

ri  (2.1c) 

a-i  a-i 

=  (Pi  +  Qi)  Uitt(2Ri)  -  (P2  +  Q2)U2  tt(2R2) 
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where  &  denotes  internal  energy  and  W  here  denotes  energy 
per  unit  mass  added  to  the  medium  from  external  sources. 

If  all  quantities  are  continuous,  the  limiting 

process  rt  -*•  ra  may  be  carried  out  in  the  usual  way  to  pro- 

2 

duce  the  differential  equations 

It  *v  ir<2R)a"1  H]  -  0 

|t^V  U  (Hi  1  =  “  Htp  +  ^1  *  7r(2R) 

It [^7 u2  +  &  -  n)1^211)  Hi =  -  (p  +  Qiui 

a-i 

y  tt(2r) 

This  set  of  equations  may  be  simplified  by 
noting  that  the  mass  equation,  Eq.2.2a,  may  be  immediately 
integrated  to  give 


(2.2a) 

(2.2b) 

(2.2c) 


-j  'v 

V  7r(2Rl  3F  “  k 

where  k  is  an  integration  constant.  At  t  =  0,  R  =  r, 
V  =  V0  so  that 

k  =  (2r)a  1  tt/Vo 

and 


1  R 

V  R 


a-i 


dR 

3F 


a-i 


(2.3a) 


Substitution  in  Eqs.2.2b  and  2.2c  then  yields 


a-i 

-  vo  |^(p  +  Q)  •  (§)  (2.3b) 
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and 


1^(1  U2  +  e)  =  -  Vo  ^[(p  +  Q)U]  (|)  +  (2.3c) 

Further  manipulation  then  yields  the  familiar 
equations  of  motion 

V-Vo^cf)0"1  (2.4a) 

If  -  "  v  Ir<p  +  Q)  (2-4b) 

f|--(P  +  Q)f|  +  |£  (2.4c) 

Together  with  the  definition 

-  U  (2.4d) 

an  equation  of  state  describing  material  behavior 


P«f(V,s)  (2.4e) 

and  an  equation  defining  the  dissipative  term  Q  in  terms 
of  the  other  variables,  these  are  then  sufficient  to  de¬ 
fine  a  solution  subject  to  given  initial  and  boundary  con¬ 
ditions. 

The  differential  equations,  Eq.2.4,  may  be  dif¬ 
ferenced  in  a  variety  of  ways.  It  is  of  course  desirable 
that  the  difference  equations  correspond  to  the  differential 
equations,  Eq.2.4,  to  second  order,  i.e.,  a  Taylor  series 
expansion  of  the  difference  expression  about  a  point  should 
correspond  to  the  differential  expression  in  terms  of  first 
order.  However,  it  is  also  considered  desirable  that  the 
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difference  equations  should  provide  second  order  analogs  of 
Eqs.2.2  and  2.3,  i.e.,  explicitly  conserve  mass,  momentum, 
and  energy,  also  to  second  order.  This  restricts  the  choice 
of  difference  equations  somewhat. 

The  continuum  is  divided  into  finite  segments  Ar  in 
(Lagrangian)  space,  and  At  in  time.  The  usual  notation 


Xn  =  X(mAr,  nAt) 
m 

is  used.  Providing  all  quantities  are  smooth  and  continu¬ 
ous,  and  vary  sufficiently  slowly  over  the  finite-difference 
mesh  (i.e.,  the  mesh  is  chosen  small  enough),  quantities 
may  be  averaged  over  each  mesh. 

Considering  the  segment  mAr  to  (m  +  l)Ar,  Eq. 

2 . la  thus  becomes 


d 


a-i 

tt2 

aV 

m+i  ya 


or  using  central  time  differencing 


0 


ac-i 


tt2 


nil 


aV 


mi  l  y* 


nil  a'  nil  a 

a-i 

tt2 

a  a 

1 

H 

3 

%e 

i 

H 

=  i 

Pi 

- 

aV  1 

mil  /z 

nfi  /2 
At  7 


0 


which  is  a  second  order  analog  of  Eq.2.2a. 

Neglecting  second  order  terms,  this  may  also  be 

written 


yri-i 
mil  y2 


n+i  a-i 

) 


(R, 


n+i 

nH-i 


R1"*1 ) 

m 


V° 

mil  y2 


a-i 


(r, 


mil 


rJ  (2-5a) 


This  is  a  second  order  analog  of  Eq.2.3a.  Equation  2.5a 
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could  have  been  written  down  directly,  since  mass  is  con¬ 
served  by  definition  of  the  Lagrangian  mesh.  However, 
proceeding  in  the  same  way,  equations  which  are  second 
order  analogs  of  Eqs.2.1b  and  2.1c,  2.2b  and  2.2c;  and  2.3b 
and  2.3c  may  be  derived.  The  analogs  to  Eq.2.3b  and  2.3c 
are 


n+i 


n 


Unrt-1/2  —  Unr+i  /  2 
_  n+1/2 


VgrH/a 


At 


n+1/2  V*1/Z\ltf+/Z 


a- 


<rnri-i  ~  rJ  L 

n+l/2  n+l/2' 


m+i 


+  Q 


m+i 


rm+i/2 


n+i/2 


a-i  _ 


-  P_ 


m 


<4 


m 


rm+i/2 


(2.5b) 


T  |  n+i  \2  1  I  \  2  .  n  n+1/2 

J  I'Wi/aJ  -  J  (uS+i/s)  +  eJSi/.,  -  AW®h/* 

.n+1/2 


At 


Vm+i  /  2 

<rnrt-i  "  rJ 


I  n+i/2  n+1/2  \  n+i/2 

Pnr+i  +  ^nr+i  j  Um+i 


R 


n+1/2 

m+i 


a-i 


n+1/2  n+1/2 

P  +  Q 

m  m 


^11^1/2 
lri-x/2  /Rtri'1/S  \“'1  I 


(2.5c) 


u 


m 


nr+i/2 

Neglecting  only  second  order  quantities,  manipula- 


tion  of  these  equations  yields 

n+i  n+i 

n+i  . 

n+i  _ 

V  =  v 0 

fRnr+i  -  Rm  | 

Rtrr+i/2 

m+i  /  2  m+i  /  2 

rm+x  -  rm  J 

rnrl-i/2  j 

a-i 


(2.6a) 
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rri-i  xi 

UnrH/a  ~  Unrfi/a 

At1*1/2 


V 


ttfl/2[PttH-i 
m+i  /2 


rrt-i  /a 


pr/a  + 


n+i  ya_  _n+i  /a 


Rnrfi 


(2.6b) 


piH-i  on 

nrt-1/2  -  mfi/2 


n+i/2  ,  .  n+i/21 
pnH-i  /2  ^  Qvrrf  1  /2 


+  AW, 


rH-i/2 

wH/2 


V, 


.rt+i 


m+i  /s 


It  may  be  easily  verified  that  these  correspond 
to  Eqs.2.4  to  second  order. 

If  quantities  vary  sufficiently  slowly  over  the 
mesh,  then,  to  second  order 


Xmfi/2  +  Xm^ 


The  equations  may  be  rendered  explicit  by  moving 
Eq.2.6b  back  half  a  mesh  in  space  and  time,  and  expressing 
Eq.2.6a  in  its  equivalent  form 


n+i 

Vm+i  /a 


n-fi 


n+i. 


=  V° 


”rt'1/a  (U  )“  -  (r_)“ 

v  nri-i  '  ni 


(2.7a) 


n+i  /2  n-i  /2 

Urn  —  Um  _ Pm-i/a  ~  pmfi /2  +  Qm-i/e  ~  Qrcri-i /2 _ 

(2.7b) 
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X 


rri-i 

^m+i  /2  ~  ® 


n 


rt+i 


n 


rrt-i 

itH-i/2  =  i(Pm+i/2  +  Pm+i/2  +  ^nH-i/2  +  ^ntfi/21 


(2.7c) 


f  vn  —  v  n^*1 

trrfi  /2  m+i  ^/2 


trl-i/2 

j  +AWnrt-i/2 


These  equations  are  then  taken  together  with  the 
difference  forms  of  Eqs.2.4d  and  2.4e 


n+i 


R 


m 


-■£ 


At 


rrt-i  Jl 


U, 


rrt-i  /z 


m 


(2 • 7d) 


P1^1 

tid-i  ^/2 


nfi 

f  (  V 

rrrt-i  /2 


n+i 

e  ) 

m+-i/2 


(2 • 7e) 


3 

Following  von  Neumann  and  Richtmyer,  a  smooth 
solution  is  ensured  by  a  quadratic  artificial  viscosity 
term.  Provision  is  also  made  for  a  linear  artificial  vis 
cosity  term. 


n  - ;p -  I  T 2  |AU|  +  c  c 

Vnrt-v=  1 


m+V 

where 


n-i 

m+i /s 


AU 


n-i  /z  n-i  /z 


AU“UmH  "Um 


and 


c  2  =  dp/dp 


(2.7f) 


represents  sonic  velocity.  Constants  q  and  C  are  normally 
set  to  zero  when  AU  <  0. 

The  conditional  stability  of  the  explicit  equa¬ 
tions  is  well  known.  An  approximate  stability  criterion. 
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for  nearly  isentropic  processes,  derived  in  a  manner  analo- 

3 

gous  to  that  of  von  Neumann  and  Richtmyer  has  proved  use¬ 
ful  in  many  cases. 


rri-i  /2 
At 

nrt-i  /a 


R 


n 

m+i 


cJVi/2  (1  +  20  +  ^2|AU| 


<2. 7g) 


the  minimum  over  all  m  being  used  to  advance  the  calcula¬ 
tion. 


Equations  2.7f  and  2.7g  are  not  correctly  centered, 
but  this  should  not  entail  serious  consequences  since  all 
quantities  are  necessarily  varying  slowly.  (In  the  deriva¬ 
tion  of  Eq.2.7g,  quantities  are  taken  as  locally  constant.) 
Correct  centering  of  these  quantities  in  Eqs.2.7b  and  2.7c 
is  much  more  important . 

Equations  2.7  are  similar  to  those  used  by 
Grandey, ^  except  for  the  centering  of  Q.  It  may  also  be 
noted  that  the  equations  reduce  identically  to  those  of 
Richtmyer'*  with  correct  Q  centering,  when  equal  space  and 
time  increments  and  uniform  initial  specific  volume  are 
used. 

If  all  quantities  are  known  at  time  n  (or  initi¬ 
ally  at  time  zero,  corresponding  to  initial  conditions) 
Eq.2.7g  is  first  used  to  determine  the  maximum  allowable 
time  interval  to  advance  the  calculation,  followed  in  turn 
by  Eqs .2. 7b,  2.7d,  2.7a,  and  2.7f  to  obtain  u”*1/2,  R1*1  , 

V 0+1  and  Q1^1  ,  respectively  over  all  m.  Equations  2.7c 

and  2.7e  are  then  solved  as  simultaneous  equations  for 
rri~i 

P  for  all  m.  This  process  is  then  repeated  stepwise  to 
produce  the  desired  solution. 

Boundary  conditions  are  introduced  through  the 
initial  conditions,  e.g.,  layers  of  different  materials  are 
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introduced  by  specifying  appropriate  initial  specific  volumes, 
and  providing  for  different  equations  of  state  in  the  appro¬ 
priate  meshes.  Free  surfaces  may  be  introduced  by  setting 
the  initial  density  and  pressure  to  zero  in  the  mesh  adjacent 
to  the  surface.  In  practice  it  is  somewhat  more  convenient 
to  set  the  density  and  pressure  in  the  mesh  outside  the  bounda¬ 
ry  to  zero  at  each  time  cycle  by  using  (for  a  right-hand 
boundary) 


cv*  - 

At1*1/2  +  £tn-1/2 


n 


=  V 


n 


M-i/s 


M-i/s  + 

Rn  _  Rn 

rM  KH-i 


(2.8) 


in  place  of  Eq.2.7b.  This  saves  carrying  an  extra  mesh  out¬ 
side  each  free  boundary. 

In  addition,  a  "spall  and  join"  routine  is  pro¬ 
vided,  which  operates  as  follows.  When  the  stress  reaches 
some  predetermined  critical  negative  (tensile)  value,  the 
material  is  allowed  to  part  and  form  two  new  free  surfaces. 
If  at  some  later  time  these  free  surfaces  attempt  to  cross, 
then  the  correct  interface  condition  is  again  used.  At  an 
unbonded  interface,  of  course,  the  critical  stress  is  zero, 
while  in  the  interior  of  a  material,  the  critical  stress  is 
the  fracture  stress. 
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SECTION  III 


CHARACTERISTICS  METHOD 


2 

The  method  of  characteristics  is  well  known. 
Restricting  the  discussion  to  one- dimensional  slab  symmetry 
and  isentropic  processes  (Q  -  0),  the  Lagrangian  mass  and 
momentum  conservation  equations,  Eqs.  2.4a  and  2.4b  re¬ 
duce  simply  to 


dV  „  dU 
=  V°  3r 


(3.1a) 


dU 

3E 


(3.1b) 


while  the  energy  equation,  Eq.2.4c,  becomes  simply 


dS 

3t  = 


0 


(3.1c) 


where  S  is  the  specific  entropy.  Using  Eq. 3.1c, the  equa 
tion  of  state,  Eq.2.4e  reduces  simply  to 


P  =  f(V) 


(3. Id) 


Providing  that  the  Jacobian 


we  can  interchange  dependent  and  independent  variables  in 
Eqs. 3. la  and  3. lb: 


dr 

dU 


(3.2a) 


dr 

dP 


(3.2b) 
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Introducing  new  functions 


c 


2 


L 


Eqs.3.2  transform  to 


dr 


c 


dt 

Sl 


(3.3a) 

(3.3b) 


(3.4a) 


dr 


(3.4b) 


The  system,  Eqs.3.4,  has  fixed  characteristics  in  the  U,L 
plane  given  by  the  family  of  curves 


U  +  L  «  Constant  (3.5) 

Writing 

U  —  L  =  a 
U  +  L  =  3 

Eqs.3.4  become 

^  +  c  |§  "  0  (3.6a) 

(3.6b) 

and  the  characteristic  directions  in  the  r,t  plane  are 
given  by 

|f-±c  (3.7) 


Here  c  may  be  identified  as  the  Lagrangian  sound  speed. 

While  analytical  solutions  can  be, easily  ob¬ 
tained  for  simple  wave  systems  (e.g..  Ref. 21),  it  is  more 
convenient  to  use  a  simple  numerical  scheme  for  problems 
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involving  more  complex  wave  interactions  and  reflections. 

The  solution  is  represented  by  a  finite  number  of 
characteristic  lines  defined  by  Eqs.3.5  and  3.7.  Examples 
are  shown  in  Figs. 3.1  and  3.2  in  which  the  numbering  con¬ 
vention  which  has  been  used  is  illustrated. 

Positive  characteristics  I+  and  negative  charac¬ 
teristics  I  are  numbered  consecutively.  Intersection 
points  of  characteristics  I*  and  I~  in  the  r,  t  plane  are 
denoted  Sm  ,  and  quantities  at  each  point  Sm  n  are  similarly 
subscripted  (  )m 

If  quantities  are  known  at  a  previous  point  Sm_1  , 

S  „  Eqs.3.5  and  3.7  allow  quantities  at  S  „  to  be  evalu- 

m,n-i'  ^  ^  m,n 

ated,  i.e., 

<U  +  L)m,n-i  -  +  L>m,n  <3'8a> 

-  <U-LV,n  <3'8b> 

P.V  and  c  at  S  _  may  then  be  found  from  L  „  by  Eq.3.3  once 
*  m  j  n  trij  n 

a  suitable  equation  of  state,  Eq.3.1d,  has  been  inserted. 

Then 


r 

m,n 

t 

m,n 


‘m.n-1 


1  , 
T.  <° 


+  C  ) 

m,n  m,n-i' 


(3.9a) 


r  —  r 
m,n  m-i ,  n 

t  -  t 
m,n  m-i,n 


\  (c 


m,n  +  cm-i ,n^ 


(3.9b) 


In  this  way  a  complete  solution  U(r,t),  P(r,t)  and  V(r,t) 
may  be  constructed  from  the  initial  conditions  defining  these 
quantities  at  t  =0. 


Boundaries  require  special  treatment.  At  a  free 
surface,  P  =  0  and  thus  c  =  cQ,  L  =  L0  and  V  =  V0  where  (  )o 
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refers  to  conditions  at  t  -  0  at  the  boundary.  Thus  Unim 
may  be  found  from  Eqs.3.8a  or  3.8b  directly  for  a  right- 
hand  or  left-hand  boundary,  respectively.  In  Lagrangian 
coordinates,  the  position  of  the  boundary  is  fixed  and 
r  >=  r  „  .  Thus,  t  may  be  found  from  Eqs.3.9a  or  3.9b 
for  a  right-hand  or  left-hand  boundary,  respectively. 

At  an  interface  between  different  materials,  the 
conditions  are  more  complex.  Denoting  quantities  in  the 
second  material  (considered  to  lie  to  the  right  as  in  Fig. 3.1) 
by  (')',  we  have  the  condition  that 


i 

P  ■  P 
m,m  m,m 

By  virtue  of  Eqs.3.1d  and  3.3  we  may  write 


(3.10) 


P  -  g(L) 

so  Chat  the  condition  becomes 

8  (L )  _  „  -  g'(I.')_ 

nijiri  'THjtn 

Furthermore  one  has  the  condition 

U  *  u' 
m,m  m,m 

In  the  left-hand  material,  Eq.3.8a  applies,  i.e.. 


(3.11) 


(3.12a) 


(3.12b) 


(U  +  L)  =  (U  +  L)  ,  _  (3.12c) 

'  nijin  7  m-i  jin  x 

In  the  present  example,  a  simple  wave  is  propagating  to  the 
right  into  the  second  material  (Fig. 3.1).  Thus,  (U  —  L)  is 
constant  everywhere,  and 

<«’  -  <U'  (3-12d) 

Equations  3.12  are  four  simultaneous  equations  for 
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the  four  unknowns  L  L*  ,  U  ,  and  u'  .  Similar  equa- 

lOjin  m^m  m,vc\ 

tions  may  be  easily  written  for  the  case  of  a  wave  moving  to 
the  left  being  refracted  at  an  interface. 


If  surface  positions  are  desired,  the  simplest 
scheme  is  to  accumulate  displacements  at  each  mesh  point. 


i.e. 


R_ 


mm 


mm 


rn 

I 

k=i 


1 

1 


(U 


kk 


Uk-i,k-i^ 


^kk  “ 


tk-i  ,k-i 


) 

(3.13) 


The  calculation  procedure  can  be  mechanized  to  al¬ 
low  use  of  a  digital  computer.  The  most  difficult  task  is  to 
write  the  logic  necessary  to  determine  the  configuration  of 
the  characteristic  mesh  in  the  r,t  plane.  Providing  that 
the  form  of  the  characteristic  mesh  is  known  beforehand,  this 
can  be  easily  done.  Two  approaches  have  therefore  been  used. 
Programs  have  been  written  (for  the  Royal  McBee  LGP  30 
computer),  covering  simple  interactions  as,  for  example,  re¬ 
flection  of  a  wave  fan  at  a  free  surface,  or  crossing  of  two 
wave  fans.  Also,  Eqs.3.8  through  3.12  have  been  programmed 
so  that  the  computer  may  be  operated  in  a  semi -manual  manner, 
i.e.,  quantities  at  points  Sm^n_1  and  Sm_1^rtare  entered 
manually  at  the  keyboard  of  the  machine,  and  quantities  at 
Sm  n  are  then  computed  automatically  upon  selection  of  the 
correct  equations,  e.g.,  interior  point,  free  surface,  etc. 

In  this  way  solutions  may  be  built  up  quite  rapidly  with  a 
minimum  of  complex  computer  logic. 

The  foregoing  has  been  restricted  to  isentropic 
processes,  and  thus  strong  shocks  cannot  be  handled  by  this 
analysis.  It  is  possible  to  extend  the  method  to  include 
entropy  variations,  in  which  case  a  third  characteristic 
direction  is  found,  representing  particle  paths,  along  which 
the  entropy  is  constant.  Strong  shocks  are  then  treated  as 
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internal  floating  boundaries  across  which  the  Rankine-Hugoniot 
relations  are  used. 

However,  for  weak  shocks  across  which  the  entropy 
change  is  negligible,  it  is  possible  to  use  the  weak  shock 
approximation  of  Courant  and  Friedrichs, ^  within  the  frame¬ 
work  of  the  foregoing  isentropic  theory. 

Thus,  when  two  positive  characteristics  meet  (as 

J. 

in  Fig. 3. 2),  say  Ik  ,  Im  ,  a  shock  wave  forms  representing 
a  jump  in  variables  f rom  (  to  (  )m  The  velocity  of 

the  shock  and  hence  its  slope  in  the  r,t  plane  is  approxi¬ 
mated  by 

||  =  i(c,  +  c  )  (3.14) 

dt  2V  k,n  m,n'  v  ' 

When  several  positive  characteristics  combine,  the 
subscripts  k,m  refer  to  conditions  on  the  I+  characteristics 
of  lowest  and  highest  number  respectively,  involved  in  the 
shock  formation.  Points  on  the  characteristic  mesh  involving 
a  shock  are  labeled  with  three  subscripts,  k,  m,  and  n,  as 
shown  in  Fig. 3. 2  for  a  right  going  shock. 

The  weak  shock  may  be  regarded  as  a  group  of  co¬ 
incident  characteristics,  and  a  point  Si  ■  —  may  therefore 

K.  j  in  j  n 

be  thought  of  as  representing (m-k)  coincident  points  n, 

J>k+i  n  •  •  •  Sm  .  Thus,  interaction  of  a  right-going  shock 
with  a  negative  characteristic  or  a  surface  may  be  handled 
(in  the  computer)  as  a  series  of  (m-k)  characteristic  inter¬ 
actions,  using  Eqs.3.8.  However,  Eq.3.14  is  substituted  for 
Eq.3.9a,  to  obtain  r  and  t. 

A  left-going  shock  may  be  treated  in  the  same  way; 
in  this  case  the  last  two  subscripts  refer  to  the  first  and 
last  i”  characteristics  involved  in  the  shock  formation.  The 
system  requires  use  of  four  index  registers  in  the  computer. 
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a  nonzero  value  in  the  first  register  signalling  the  existence 
of  a  right-going  shock,  and  a  nonzero  number  in  the  last  regis¬ 
ter  signalling  the  existence  of  a  left-going  shock. 
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SECTION  IV 


EQUATION  OF  STATE 


4.1  Introduction 

Strictly  speaking,  a  real  solid  does  not  in 
general  possess  a  simple  equation  of  state  in  the  sense 
that  there  exists  a  unique  single-valued  function  relating 
the  variables  characterizing  the  state  of  the  material.  Due 
to  effects  such  as  strain  hardening,  phase  changes  and  the 
possibility  of  nonequilibrium  states,  the  relation  among 
state  variables  is  dependent  on  the  entire  history  of  the 
material  up  to  the  time  under  consideration.  However, 
with  certain  necessarily  gross  assumptions  regarding  the 
initial  state  of  the  material  before  commencement  of  the 
process  under  analysis,  and  the  allowable  changes  in  the 
material  during  the  process,  it  is  possible  to  construct  a 
relation  among  state  variables  at  a  particular  time,  which 
is  a  function  of  the  detailed  path  of  the  process  up  to  that 
time.  The  term  constitutive  equation  may  be  applied  to 
cover  such  relationships.  Once  such  an  equation  has  been 
constructed,  it  is  quite  possible  to  use  it  in  a  numerical 
solution  of  an  initial  value  problem,  particularly  in  a 
Lagrangian  finite-difference  scheme,  since  it  is  possible 
for  the  computer  to  keep  track  of  the  detailed  history  of 
each  mass  point. 

At  very  high  pressures,  it  is  usually  assumed 

that  the  deviatoric  stress  component,  which  the  material 

may  support,  is  negligible  in  comparison  with  the  dilatational 

stress.  If  equilibrium  is  attained  very  rapidly  (about  10~7 

seconds)  and  if  wave  profiles  are  steady  in  time,  then  a 

simple  single-valued  entropic  equation  of  state  may  be  found 

7-9 

from  measurements  on  strong  shock  waves. 
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It  is  here  assumed  that  an  extrapolation  of  this  equa¬ 
tion  to  low  stresses  is  still  valid  for  the  dilatational  com¬ 
ponent.  This  assumption  is  investigated  later  in  this  report 
by  comparison  with  static  compressibility  and  thermodynamic 
data  at  zero  pressure. 

If,  further  it  is  assumed  that  the  deviatoric  com¬ 
ponent  of  stress  is  limited  by  a  yield  condition  and  that  the 
dilatational  component  of  plastic  strain  is  zero, ^  it  is 
then  possible  to  deduce  a  relation  between  the  longitudinal 

components  of  stress  and  strain  from  weak  shock  wave  measure- 

17  18 

ments.  For  a  few  aluminum  alloys  Fowles  and  Lundergan 
have  shown  that  the  resulting  relation  is  comparable  with 
static  stress-strain  relations  obtained  under  conditions  of 
uniaxial  stress  (engineering  tensile  test),  indicating  that 
equilibrium  is  also  reached  very  rapidly  in  these  materials. 

4.2  High  Pressure  Region 

The  high  pressure  equation  of  state  has  been  under 

7  8  9 

investigation  for  some  time.  '  *  The  experimental  method 
involves  inducing  a  plane  shock  wave  propagating  into  un¬ 
stressed  material  under  investigation  by  explosive  or  other 
means,  and  determining  the  shock  velocity  and  material  ve¬ 
locity  behind  the  shock.  If  the  shock  wave  profile  remains 
constant  In  time,  and  equilibrium  is  achieved  very  rapidly 
behind  the  shock,  the  Rankine*flTIgoniot  relations  then  allow 
the  compression,  stress,  and  energy  change  across  the  shock 
to  be  related  to  the  shock  velocity  Uc  and  the  material  ve- 
locity  Up  behind  the  shock  relative  to  the  material  ahead 
of  the  shock. 


(4-la) 
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pcL  '  :  jfi  r 


UsUP 


_  0$  Of 


«=> 


H 


j/vT f*tVl>y 


eH  “  eo  ‘  7  rHv,° 


P„(V0  -  V) 


(4.1b) 

(4.1c) 


Here  (  )u  refers  to  conditions  behind  the  shock, 
and  (  )0  refers  to  ambient  conditions  in  the  unstressed 
material  ahead  of  the  shock. 

The  locus  of  P„(V)  points  determined  from  many 

. — — —  H 

such  experiments  with  different  shock  strengths  is  termed 
the  shock  Hugoniot. 

For  many  materials  it  is  found  experimentally  that 

U  is  a  linear  function  of  Uc  over  a  wide  range,  expressible 
S  r 

by  the  relation 


Ug  =  c  +  s  Up 


(4.2) 


Using  Eqs.4.1  and  4.2,  PH(V)  becomes 


c2(Vo  ~ 


-  s(V0 


(4.3) 


The  shock  Hugoniot  determined  in  this  way  repre¬ 
sents  states  which  may  be  reached  by  a  shock  process  from 
the  unstressed  state  at  ambient  conditions.  It  is  de¬ 
sirable  to  find  a  relationship  between  P,V,  and  C  for 
states  neighboring  the  shock  Hugoniot,  which  may  be  reached 
by  neighboring  adiabats  or  isotherms.  Additional  informa¬ 
tion  is  required.  The  usual  procedure  is  to  assume  a  func¬ 
tional  form  for  the  equation  of  state  based  on  solid  state 
theory.  This  has  been  done  by  Walsh  et  al.^  The  theory 
will  be  briefly  reproduced  to  facilitate  later  discussion. 
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The  equation  of  state  is  assumed  to  take  the  form 

12 

of  the  Mie-Gruneisen  equation 

(P  -  PK)  -  7 

Here  (  )R  refers  to  t 
the  Gruneisen  ratio,  is  the  ratio  of  the  thermal  component 
of  pressure  to  the  lattice  vibrational  energy  density. 

The  assumption  is  now  made  that  at  very  high  pres¬ 
sures  the  deviatoric  stress  which  the  material  may  support 
is  negligible  in  comparison  with  the  dilatational  stress, 
and  the  state  of  stress  may  be  simply  characterized  by  a 
pressure.  Under  this  assumption,  Eq.4.4  may  be  written 


e  -  e 


K 


*  0  °K  isotherm,  and  7, 


<PH  -  PK> 


(4.5) 


While  attempts  have  been  made  to  estimate  7  at 

9  10 

very  large  compressions  from  experimental  data,  *  Walsh 

et  al.  use  a  relation  between  7  and  P„  derived  from  solid 

12  ^ 
state  theory.  Slater  has  given 


V  'Sv? 

7  “  “  2  3~7p'T 
3v  (PKJ 


2 

1 


13 

while  Dugdale  and  MacDonald  have  given 


(4. 6a) 


7 


v  Iv*  (pKv2/3>  _  1 
2  iv  (pkvV3)  3 


(4.6b) 


Inserting  Eq.4.6a  or  Eq.4.6b  together  with  Eq.4.1 
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into  Eq.4.5,  and  noting  that 


SK  ‘  f°  PK<V)  dv 
V 

results  in  a  differential  equation  for  P^,  which  may  be 
solved  (in  this  case  numerically)  with  suitable  initial 
values  P^,(dP/dV)^  ,  V*  and  e'  at  some  known  point  on  the 
0  °K  isotherm.  The  quantity  7(V)  may  then  be  determined 
from  Eq.4.6a  or  Eq.4.6b. 

In  order  to  obtain  suitable  initial  conditions, 
it  is  necessary  to  introduce  additional  assumptions.  At 
low  pressures  the  deviatoric  stress  which  the  material  may 
support  is  not  negligible.  Since  Eq.4.4  applies  only  to 
pure  dilatation,  shock  Hugoniot  data  obtained  at  low  pres¬ 
sure  cannot  be  used  to  proceed  from  known  initial  conditions 
at  zero  pressure.  One  might  make  some  gross  assumptions  in 
order  to  arrive  at  suitable  initial  conditions  at  suffici- 

Q 

ently  high  pressure,  or  it  is  necessary  to  make  some  as¬ 
sumptions  regarding  extrapolation  of  the  shock  Hugoniot 
data  to  zero  pressure.  The  latter  course  has  been  used 
by  Walsh  et  al.^  It  is  assumed  that  the  Ug  versus  Up 
curve  may  be  linearly  extrapolated  to  Up  =0.  For  many 
materials,  the  intercept  c  corresponds  quite  closely  to  the 
adiabatic  bulk  sound  speed  deduced  from  static  compressibility 
measurements. 

Initial  conditions  then  become  P^  =  0,(dP/dV)OK> 

V  ,  and  g  ,  where  (  )Q  refers  to  conditions  at  the 
foot  (P  =  0)  of  the  0  °K  isotherm.  It  is  possible  to  obtain 
(  e0  -  6  or)  by  integrating  the  Debye  specific  heat  curve. 

The  quantities  VOK  and  (dP/dV)oK  may  be  estimated  from  Eq. 

4.4,  i.e.  at  P  =  0,  V  VqK; 
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VoK  "  7oK 


(4.7a) 


eH^VoK^  ”  SoK 

k) 


dP 

-5V 


OK 


+  PH(V°k) 


(4.7b) 


a  eH 
3v 


oK 


} 


These  may  be  evaluated  if 
suitable  values  of  7oK  and 
(S7/SV)0^  can  be  found. 
Assuming  (c^/dV)  constant  in 
the  small  interval  VoK  to 
V0  we  have 


(I?)  -  ( 

3?) 

(4.8a) 

ldV/  OK  l 

4 

7oK 


7 

o 


(4.8b) 


Assuming  the  Slater  relation,  Eq.4.6a,  or  Dugdale- 
MacDonald  relation,  Eq.4.6b,  to  be  valid  along  an  adiabat 
though  VQ  allows  70  and  (cfy/dv)  0  to  be  determined.  Using 
Eq.4.3  and  Eq.4.6a 

70=2s-|  (4.9a) 
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(!?)0  '  ^  <■**  +  3  s>  (4-9b> 

Alternatively,  using  Eq.4.3  and  Eq.4.6b 

7o  =  2s  —  1  (4. 10a) 

(!v)0  *  ^s2  “  Is  +  ^  (4.10b) 

It  may  also  be  noted  that  differentiation  of 
Eq.4.4  at  constant  volume  yields 


■Yo 


V 


Vo  I  Sp  1  I  dV  I 

S  w)sUip 


70 


Vo*  K 


(4.11) 


where  is  the  specific  heat  at  constant  pressure,  K  is 

the  adiabatic  bulk  modulus  and  x  is  the  volumetric  thermal 

expansion  coefficient,  all  at  zero  pressure. 

7  8 

Walsh  et  al.  ,  and  McQueen  and  Marsh  have  given 

analytical  fits  to  Pjj(V)  and  to  7(V)  computed  using  the 

Dugdale -MacDonald  relation,  Eq.4.6b,  and  the  thermodynamic 

definition  of  7Q,  Eq.4.11,  for  a  large  number  of  metals. 

Equation  4.4  may  be  written  as 


(P  -  PH)  -  7 


(4.12) 


It  is  often  more  expedient  for  computational  pur¬ 
poses  to  express  Eq.4.12  in  the  form 
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P  «  M(V)  +  ft  •  N(V) 


(4.13) 


where  in  this  case 

M(V)  -  P„(V)  [l  -  9f  V2_^_V  J 
N(V)  -  ^ 

It  is  generally  convenient  to  fit  M  and  N  with  polynomials 
in  (Vo/V).  Equation  4.13  may  then  be  used  directly  as  an 
equation  of  state  at  high  pressures. 

The  accuracy  of  the  pressure  offset  from  the  ex¬ 
perimental  shock  Hugoniot  (P  —  P^)  is  equal  to  the  accuracy 
with  which  7  has  been  determined.  This  depends  in  turn  on 
the  validity  of  the  linear  extrapolation  of  the  U  versus 
Up  data  to  Up  =  0,  and  on  the  validity  of  the  relation  be¬ 
tween  7  and  P^. 

In  order  to  investigate  the  uncertainties,  Eq. 

4.5  was  reintegrated  using  both  the  Slater  and  Dugdale- 
MacDonald  relations  for  7,  with  particular  attention  to  the 
low  pressure  region.  In  addition,  since  7  was  found  to  vary 
only  slowly  with  V,  an  integration  was  performed  with  7  held 
constant  at  70  defined  by  Eq.4.11. 

Numerical  results  for  copper  are  shown  in  Table  I 

for  relatively  low  pressures.  Also  shown  are  the  relevant 

initial  conditions.  U  versus  Un  data  obtained  by  McQueen 
8  s  9  r 

and  Marsh  and  Altshuler  were  fitted  by  Eq.4.2.  For  the 
case  of  copper,  the  resultant  value  of  c  is  in  excellent 
agreement  with  the  adiabatic  sound  speed  computed  from 
Bridgman's  static  compressibility  data.  Hugoniot  pres¬ 
sures  P„  listed  in  the  table  are  computed  from  the  resultant 

7 

relation,  Eq.4.2.  Walsh  et  al.  noted  that  70computed  from 
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the  Dugdale-MacDonald  relation,  Eq.4.10a,  are  generally  in 
slightly  better  agreement  with  thermodynamic  data,  Eq.4.11, 
than  when  computed  from  the  Slater  relation,  Eq.4.9a.  This 
is  true  in  the  case  of  copper,  the  Dugdale-MacDonald  -y  be¬ 
ing  in  excellent  agreement  with  the  thermodynamic  yQ.  Initi¬ 
al  conditions  for  all  three  integrations  of  Eq.4.5  were 
therefore  computed  using  the  thermodynamic  7  ,  and  the  value 
of  (d7/dV)0  from  Eq.4.l0b. 


The  differences  in  7,  and  hence  in  the  pressure 


offset  using  the  Dugdale-MacDonald  relation  and  the  Slater 
relation  respectively,  amount  to  about  8  per  cent  at  100  kb 
and  increase  to  18  per  cent  at  900  kb.  However,  since  the 
thermal  pressure  (PH  —  PR)  over  the  range  reaches  a  maxi¬ 
mum  of  about  9  per  cent  of  P^  at  300  kb,  the  computed  values 
of  PK  agree  within  2  per  cent  up  to  300  kb.  Since  the 
Slater  7  increases  slightly  while  the  Dugdale-MacDonald  7 
decreases  over  this  range,  the  constant  7  gives  results  in¬ 
termediate  between  the  two. 


Isentropes  PA  computed  from  Eq.4.12  by  setting 


Vo 


P  dv 


V 


are  also  listed.  The  offsets  of  the  isentropes  from  the 
Hugoniot  (P^  —  P^)  reach  a  maximum  of  2  per  cent  at  250  kb, 
and  are  therefore  negligible,  a  result  which  may  be  sur¬ 
mised  from  the  fact  that  the  entropy  change  is  of  third 
order  in  the  compression,  and  hence  negligible  in  this 
range.  To  all  intents  and  purposes,  it  is  permissible  to 
neglect  entropy  changes  and  use  the  Hugoniot  as  equation  of 
state  below  about  0.25  megabars  in  copper,  where  only  adia- 
batic  wave  propagation  is  involved.  For  greater  offsets. 

Adiabatic  is  here  used  in  the  strict  sense  that  no  heat  is 
added  or  subtracted  by  external  sources. 
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such  as  isentropes  or  isotherms  through  zero  pressure  states 
at  significantly  different  temperatures,  improvement  in 
precision  is  possible  by  the  use  of  Eq.4.12,  with,  however, 
a  probable  uncertainty  in  pressure  offset  from  the  experi¬ 
mental  Hugoniot  of  up  to  20  per  cent.  In  the  present  in¬ 
stance,  taking  7  constant  in  this  range  leads  to  results  as 
good  as  any,  with  significant  simplification  to  the  equation 
of  state. 

Two  checks  on  the  linear  extrapolation  of  Us 
versus  Up  data  are  possible.  One  may  obtain  isotherms  (P^) 
through  V0,  (i.e.,  at  298  °K)  which  are  directly  comparable 
to  Bridgman’s  static  compressibility  data.  One  may  also 
evaluate  the  internal  energy  as  a  function  of  volume 
(  £pmo)  at  zero  pressure,  which  is  directly  comparable  to 
zero  pressure  specific  heat  and  thermal  expansion  data.  In 
the  latter  case,  thermodynamic  data  at  high  temperatures 
corresponds  to  increased  specific  volume.  This  corresponds 
to  extrapolation  to  negative  Up  or  tension  states. 

Table  I  shows  that  the  offset  (PH  —  PT)  is  less 
than  3  per  cent  at  30  kilobars,  the  highest  pressure  at 
which  reliable  static  data  is  available,  and  the  difference 
between  isotherms  calculated  on  the  basis  of  the  three 
methods  is  consequently  negligible.  Comparison  with  static 
data  of  Bridgman  ^  shows  agreement  within  the  combined  ex¬ 
perimental  error. 

Comparing  values  of  Ep_o  shows  that  best  agreement 
with  thermodynamic  data  is  shown  by  the  Dugdale-MacDonald  7, 
the  agreement  again  being  within  experimental  error.  This 

O 

has  also  been  noted  by  McQueen  and  Marsh.  This  agreement 
is  very  encouraging,  since  the  offset  from  the  (extrapolated) 
Hugoniot  in  this  case  is  large. 

One  might  therefore  have  some  confidence  in  the 
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linear  extrapolation  of  the  Up  versus  Ug  fit.  The  assumption 
will  later  be  made  that  the  extrapolated  isentrope  applies  to 
the  dilatational  components  of  stress  and  strain  at  low  com¬ 
pressive  and  tensile  stresses,  and  the  above  comparisons  lend 
considerable  support  to  such  an  assumption. 

Similar  results  for  aluminum  are  shown  in  Table  II. 

7  9 

u  versus  U„  data  obtained  by  Walsh  et  al.  and  Altshuler 
s  P 

were  again  fitted  with  a  linear  expression,  Eq.4.2.  In  this 
case,  agreement  of  c  with  the  adiabatic  sound  speed  from 
static  compressibility  tests  is  not  so  good.  Also,  the  Slater 
70  from  Eq.4.9a  is  in  better  agreement  with  the  thermodynamic 
70,Eq.4.11,  than  the  Dugdale -MacDonald,  70,  Eq.4.10a.  Initial 
conditions  were  nevertheless  computed  from  the  thermodynamic 
70  and  (^7/dV) 0  taken  from  the  Dugdale -MacDonald  relation, 
Eq.4. 10b. 

Differences  in  7,  and  offsets  of  the  various  curves 
from  the  Hugoniot  are  of  the  same  order  as  for  copper,  and  the 
same  remarks  apply.  However,  comparison  of  the  computed 
298  °K  isotherm  P^,  with  Bridgman's  data  reveals  differences 
which  are  significant,  in  line  with  the  disagreement  in  c. 

Also  significant  is  the  disagreement  in  computed  £p_0  and 
thermodynamic  data.  Thus,  for  aluminum  the  linear  extrapola¬ 
tion  of  the  Ug  versus  Up  data  to  Up_Q  does  not  appear  to  be 
as  satisfactory,  and  the  description  of  the  low  pressure 
and  tensile  region  may  also  be  considered  to  be  less  satis¬ 
factory  than  in  the  case  of  copper. 

An  examination  of  the  experimental  data  indicates 

that  a  quadratic  U  versus  UD  fit,  constrained  to  pass 

s  r* 

through  the  Bridgman  sound  speed  at  Up  «  0  may  prove  more 
suitable.  By  carrying  out  a  least  squares  fit,  it  was  found 
that  the  expression 

Ug  -  c  +  s  Up  +  h  Up* 
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fitted  the  experimental  data  as  well  as  the  previous  linear 
fit,  with  the  advantage  of  agreeing  with  the  Bridgman  sound 
speed  at  Up  =  0.  Table  III  shows  resultant  data  at  low  pres¬ 
sures  calculated  with  constant  7.  In  this  case  the  computed 
PT  and  Sp_0  values  show  improved  agreement  with  static 
compressibility  and  thermodynamic  data. 

At  pressures  much  higher  than  0.25  megabars,  the 

offsets  between  the  Hugoniot  and  adjacent  adiabats  become 

larger.  Table  IV  shows  Hugoniot  data  for  copper,  computed 

from  Eq.4.2  fitted  to  data  of  McQueen  and  Marsh,  and 

9 

Altshuler  et  al.  Also  shown  are  0  ©K  isotherms,  Gruneisen 

8 

ratios  and  temperatures  computed  by  McQueen  and  Marsh  using 

9 

the  Dugdale-MacDonald  relation,  and  by  Altshuler  et  al., 
who  included  the  contribution  of  the  electron  potential. 

The  two  results  compare  very  well,  so  that  it  may  be  surmised 
that  the  effect  of  the  electron  potential  is  negligible  up  to 
1.5  megabars.  The  differences  resulting  from  the  use  of  the 
Slater  relation.  Eq.4.6a,  instead  of  the  Dugdale-MacDonald 
relation  are  also  small,  the  difference  in  7  not  exceeding 
20  per  cent  over  the  entire  range. 

4,3  Low  Stress  Region 

Shock  Hugoniots  have  been  determined  for  a  few 
materials  in  the  low  stress  region  using  the  same  experi¬ 
mental  techniques  as  in  the  high  pressure  region,  i.e.,  by 
inducing  a  plane  loading  wave  in  the  sample  under  study  by 
explosive  means  or  plate  impact,  and  measuring  wave  and 
material  velocities.  Providing  that  the  wave  profiles  are 
steady,  it  is  usually  assumed  that  the  Rankine-Hugoniot  re¬ 
lations  may  be  used  to  deduce  stress  and  specific  volume  be¬ 
hind  the  wave  from  the  wave  velocity  and  material  velocity 
behind  the  wave.  However,  in  this  case,  the  stress  so  de¬ 
duced  is  the  stress  along  the  axis  of  motion.  Since  the 
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material  may  be  expected  to  support  a  shear  stress  which  is 
not  negligible  in  this  region,  the  stresses  along  principal 
axes  perpendicular  to  the  motion  will  differ  from  the  longi¬ 
tudinal  principal  stress.  For  this  reason,  only  one¬ 
dimensional  slab  symmetrical  motion  will  be  discussed  in  con¬ 
nection  with  elastic-plastic  wave  propagation,  since  for 
cylindrical  or  spherical  symmetry,  additional  terms  appear 
in  the  equations  of  motion,  Eq.2.4. 

At  low  stresses,  two  loading  waves  are  observed, 
one  traveling  at  elastic  longitudinal  wave  velocity,  and  a 
second  slower  wave  identified  as  a  plastic  wave.  The 
second  wave  travels  into  a  stressed  and  moving  medium,  so 
that  the  mechanical  Rankine-Hugoniot  relations  become 

UP  “  Up, 

Vi  -  V  -  Vx  /  _  (4.14a) 

us  uPi 


°H  ”  ai 


<UP  -  UPx}  <Us  "  V 

Vx 


(4.14b) 


where  (  ) xrefers  to  conditions  ahead  of  the  wave  and  a 
here  refers  to  the  stress  component  in  the  direction  of 
motion.  The  locus  aH(V)  is  again  termed  the  shock  Hugoniot. 

A  number  of  workers  have  given  a  simple  theory  ot 
uniaxial  strain  for  an  isotropic  material,  based  on  linearized 
elastic-plastic  theory . ^  This  may  be  summarized  as 
follows . 


In  uniaxial  strain,  the  principal  components  of 

0,  0) ;  the  principal  components  of 

Qy>  Qy) •  It  follows  that  any  component  of 

the  deviator ic  stress  may  be  written  a (a  —  a  )  where  a  is 

x  y 

an  appropriate  constant.  It  is  assumed  that  yielding  is  not 
influenced  by  the  dilatational  component  of  stress,  but 


natural  strain  are  (e  , 

A 

stress  are  (ax. 
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1 1 

: 

v  :  0  5/7V  /  r  Ofr, 

/  u  r 

depends  only  on  the  deviator ic  components,  and  consequently 

on  the  value  of  (a  —  a  ).  The  yield  condition  may  there- 

x  y 

fore  be  written 

<°x  "  V2  =  Y  2  (4*15> 

Hooke's  Law  will  be  written 

d€i  =  da^  +  da  (4.16) 

relating  the  principal  components  of  natural  strain  in¬ 
crement  to  the  principal  components  of  the  deviatoric 
stress  increment  dc|  and  the  dilatational  stress  increment 
do  where  a  =  l/3(ax  +  2ay) .  Here  G  is  the  shear  modulus 
and  K  the  bulk  modulus. 

If  the  strain  increment  is  divided  into  elastic 
and  plastic  components. 


dei  =  de®  +  deP  (4.17) 

and  if  the  usual  assumption  is  made  that  no  volume  changes 
result  from  the  plastic  components  of  strain 

deP  +  2deP  =  0  (4.18) 

x  y  v  ' 

the  principal  components  of  the  stress  increment  and  strain 
increment  in  the  elastic  and  plastic  ranges  may  be  written 
as  in  Table  V,  which  shows  also  their  dilatational  and 
deviatoric  components. 

If  the  values  from  Table  V  are  introduced  into 
Hooke's  Law,  Eq.4.16,  the  following  relations  result: 
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Elastic  dox  =  (K  +  4/3GJ  dex  (4.19) 

Plastic  de®  -  dcx  +  ~K-g^g^G  dY  (4.20a) 

<  -  k do* ?  dY  <4-20b> 

and  on  addition  of  Eq.4.20a  and  4.20b, 

da  ■  K  d€  +  2/3  dY  (4.21) 

X  A 

The  upper  sign  in  the  above  equations  is  appropri¬ 
ate  when  do  is  positive,  the  lower  when  do  is  negative. 

X  x 

Now  for  loading  below  the  elastic  limit,  Eq.4.19 
may  be  immediately  integrated  to  give 

o  =  (K  +  4/3G)  e  (4.22) 

X  A 

where  it  is  assumed  that  the  moduli  are  constant  over  the 
elastic  range.  Equation  4.22  above  may  also  be  written 

ox  =  F  ex  (4.23) 

where  F  =  (K  +  4/3G)  may  be  termed  the  longitudinal  elastic 
modulus. 

This  relationship  holds  up  to  the  yield  point, 
Eq.4.15,  where 


XI  ■  h  Yl 

(4.24a) 

Xl  ■  25  Yl 

(4.24b) 

(  ) x  referring  to  initial  yield. 
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Beyond  this  point  the  material  behaves  according 
to  Eq.4.21.  Howevei;  it  is  not  possible  to  integrate  this 
immediately  since  the  elastic  moduli  may  be  expected  to 
vary  significantly  with  compression  in  the  plastic  range. 

The  difficulty  can  be  overcome  by  relating  the  plastic 
loading  behavior  to  the  pressure-volume  relation  for  pure 
dilatation.  In  the  case  of  pure  dilatation,  if  dv  is  the 
volumetric  strain  increment  consequent  upon  a  change  of 
hydrostatic  pressure  dP  then 

dP  =  K  dv  (A. 25) 

Notice  that  in  uniaxial  strain 

v  =  ex  (4.26) 

Suppose  that  at  the  yield  point  in  uniaxial  strain 
the  volumetric  strain  is  vx  so  that  vx  =  e  and  suppose 
that  the  hydrostatic  pressure  which  would  give  rise  to  the 
same  volumetric  strain  in  the  case  of  pure  dilatation  is 
Pj  .  Then  from  Eqs.4.24,  4.25,  and  4.26  it  follows  that 

aXi  -  PL  =  2/3  V,  (4.27) 

Now  using  Eqs .4.25  and  4.26,  if  it  is  assumed  that 
K  has  the  same  value  in  Eq.4.21  and  Eq.4.25,  Eq.4.21  may  be 
rewritten 


dcrx  =  dP  +  2/3  dY 

If  there  is  no  strain  hardening,  so  that  dY  =  0, 
integration  from  a  lower  limit  at  the  yield  point  gives 

°x  “  *Xl  *  P  -  pi 
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or,  from  Eq.4.27 


(4.28a) 


a  «=  p  +  2/3  Yi 

A 

where  a  and  P  are  values  corresponding  to  the  same  volu- 
metric  strain. 

If  strain  hardening  is  included,  the  same  integra¬ 
tion  gives 


ax  -  axi  =  F  -  PL  +  2/3  (y  -  Yi ) 
so  that  Eq.4.28a  may  now  be  written 

cx  -  P  +  2/3  Y  (4.28b) 

A 

where  Y  is  the  value  of  the  yield  function  (°x  “  <*y)  in 
uniaxial  strain  at  the  stress  o  . 

A 

Fowles^  has  further  extended  the  argument  in 
order  to  relate  Y  in  Eq.4.28b  to  the  yield  stress  in  an 
engineering  tensile  test  (uniaxial  stress)  in  the  presence 
of  strain  hardening.  However,  within  the  region  considered 
here,  the  plastic  strain  under  conditions  of  uniaxial  strain 
is  quite  small  in  comparison  with  that  possible  in  a  tensile 
test  up  to  fracture.  One  might  therefore  expect  that  strain 
hardening  effects  are  very  small  for  a  material  such  as  that 
investigated  by  Fowles,  as  is  shown  in  fact  by  his  data. 

Upon  unloading  from  a  plastic  state  (a  ,  €  ), 

At  Ac 

the  material  behaves  elastically  according  to  Eq.4.19. 

Again  assuming  K  and  G  constant  over  the  elastic  unloading 
range,  but  not  necessarily  the  same  as  at  initial  loading, 
we  have 

<ax2  -  ax>  <4*29> 
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If  there  is  no  strain  hardening  and  no  Bauschinger 
effect,  the  reverse  yield  point  may  be  reached  when 


(4.30) 


Further  unloading  will  be  plastic,  defined  by 


ox  =  P  —  2/3  Yi  (A .31) 

as  for  loading.  The  resultant  stress-strain  path  is 
sketched  in  Fig. 4.1. 

Taking  the  available  data  on  low  stress  Hugoniots, 

17 

i.e.,  that  of  Fowles  for  hardened  and  annealed  2024 

18 

aluminum  alloy,  and  that  of  Lundergan  for  a  hardened  6061 
aluminum  alloy,  one  may  plot  (ox  —  P)  versus  specific 
volume,  where  P  is  taken  from  the  isentrope  computed  in  the 
previous  section.  Fig. 4. 2  a  and  b.  Drawn  horizontal  lines  cor¬ 
respond  to  two  thirds  of  the  yield  stress  at  .2  per  cent  off¬ 
set  from  static  tensile  test  data  reported  for  each  material, 
while  the  horizontal  dotted  lines  represent  two  thirds  of 
the  yield  stress  obtained  by  fairing  a  horizontal  line 
through  the  initial  portion  of  the  engineering  stress-strain 
curves  obtained  in  a  tensile  test  for  each  material.  It 
may  be  seen  that  the  latter  provides  a  quite  reasonable  ap¬ 
proximation  to  the  Hugoniot  data  in  each  case,  within  the 
experimental  uncertainty,  although  the  experimental  points 
do  suggest  some  strain  hardening  for  the  2024  alloy.  One 
might  thus  conclude  that  Eq.4.28a  provides  a  reasonable  ap¬ 
proximation  for  these  materials  for  the  plastic  loading 
path. 

The  previous  theory  rests  on  infinitesimal  strain 
19 

theory.  Backman  has  indicated  the  order  of  error  to  be 
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expected  from  neglect  of  the  second  order  terms  in  uniaxial 
strain.  Actually,  the  foregoing  theory  does  accommodate 
some  correction,  in  that  variation  in  elastic  moduli  with 
compression  is  permitted,  the  variation  being  taken  from 
experimental  data  on  a  pure  dilatational  compression. 

Thus,  for  the  initial  elastic  loading  portion,  it 

is  assumed  that  the  strain  is  infinitesimal  up  to  the  yield 

1  ft 

point  1  (Fig.4.1).  For  the  hardest  aluminum  alloy,  the 
volumetric  compression  1  —  V/VQ  ■  0.006  so  that  this  assump¬ 
tion  seems  reasonable. 

It  is  possible  to  compare  the  elastic  modulus  F 

18 

computed  from  elastic  wave  speeds  measured  by  Lundergan 
to  F  computed  from  K  given  by  Bridgman's  static  compres¬ 
sibility  data,  corrected  to  adiabatic  conditions,  and 
Poisson's  ratio  t  measured  in  a  static  tension  test. 

The  experimental  evidence  (Fig. 4. 2)  on  plastic 
loading  suggests  that  no  great  error  accrues  from  the  use 
of  Eq.4.28a  either,  rather,  the  correction  to  the  last 
term  of  Eq.4.28a  arising  from  the  correct  treatment  according 
to  second  order  theory  is  less  than  the  sensitivity  of  the 
quoted  experiments.  Thus,  while  it  is  desirable  to  use  a 

19 

second  order  theory  throughout,  as  has  been  done  by  Backman, 
for  our  present  engineering  purposes  Eq.4.28a  appears  to  suf¬ 
fice. 

Upon  elastic  unloading  from  a  plastic  state,  use 
of  Eq.4.29  with  a  suitable  value  of  F  rests  on  the  assump¬ 
tion  that  the  change  in  strain  between  points  2  and  3  (Fig. 
4.1)  on  the  elastic  unloading  portion  is  infinitesimal. 

20 

This  seems  reasonable.  With  such  an  assumption,  Birch, 
using  second  order  strain  theory,  has  shown  that  Eq.4.29  is 
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in  fact  correct,  and  that  the  relations  among  local  elastic 
moduli  are  the  same  as  in  first  order  theory.  He  has  further 
given  the  variation  in  elastic  moduli  with  compression  as 

F  -  (1  —  2e)5/2  [*  +  2 u  -  €(11X  +  10n)]  (4.32a) 

K  *  (1  —  2e)5/2  (A  +  |ii)  (1  -  7c)  (4.32b) 

where  Vand  p.  are  the  Lame  *  constants.  Evaluating  the  Lame/ 
constants  from  static  isothermal  infinitesimal  strain  data, 

K  has  been  computed.  Corresponding  values  of  P  have  been 
found  by  integrating  K  -  —  V  |dp/dvj,  and  are  listed  in  Table 
III.  The  encouraging  agreement  with  computed  isotherms  PT 
suggests  that  Eq.4.32a  may  also  be  a  realistic  description 
of  the  variation  of  the  longitudinal  modulus  with  compres¬ 
sion,  (when  corrected  to  adiabatic  conditions) . 

No  direct  experimental  evidence  is  currently 
available  to  check  the  assumptions  inherent  in  Eq.4.31  for 
the  unloading  plastic  path,  and  caution  is  hence  necessary 
in  its  use. 

The  foregoing  theory  has  also  neglected  strain 
rate  effects.  The  agreement  of  the  quoted  data  with  the 
rate- independent  theory  given  above  indicates  that  this 
assumption  may  suffice  for  present  purposes  for  aluminum. 
However,  this  may  not  hold  true  for  other  materials.  In 
particular,  rate  dependent  behavior  may  be  expected  in 
steels  exhibiting  an  upper  and  lower  yield  point  in  simple 
tension. 

Analysis  of  wave  propagation  including  a  relaxa- 

23  24 

tion  function  leads  to  unsteady  loading  wave  profiles.  * 
Conversely,  one  might  expect  that  rate  effects  are  absent 
when  steady  loading  wave  profiles  are  observed.  Under 
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unsteady  conditions,  the  Rankine-Hugoniot  relations,  Eq. 
4.14,  can  no  longer  be  used  directly  to  obtain  stress- 
strain  information  from  measured  wave  and  rear  surface  ve¬ 
locities.  More  complex  methods  of  deducing  constitutive 
equation  data  from  measurable  quantities  must  be  evolved  to 
deal  with  such  materials.  Further  discussion  of  this  point 
is  considered  outside  the  scope  of  this  report . 


4.4  Medium  Stress  Region 

It  is  useful  to  fit  the  dynamic  adiabat  by  a 
simple  expression,  either  for  use  in  the  elastic-plastic 
constitutive  equation,  or  for  use  alone  in  approximate 
calculations  when  the  expected  stresses  are  such  that  en¬ 
tropy  changes  across  shocks  are  small,  but  elastic  behavior 
is  insignificant.  A  suitable  expression  is  the  Murnaghan 
expression  — — 

P  -  A-[  (|a)5  _  ij.  (4.33) 

The  corresponding  Bulk  Modulus  is 


K 


-  V 


dP 


(4.34) 


Since  entropy  changes  are  neglected,  we  may  use 
the  Hugoniot  data  directly  in  fitting  the  constants  in  Eq. 
4.33.  For  the  linear  Ug  vs  Up  fit,  the  Hugoniot  P(V)  is 
given  by  Eq.4.3.  The  corresponding  Bulk  Modulus  is 


(4.35) 


or 
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K  -  -  c*  f(^) 

A  simple  means  of  fitting  the  constants  of  Eq. 

4.33  to  Eq.4.35  suggests  itself.  From  Eqs.4.34  and  4.35 

-  i  log  (~)  -  k  -  log  f(^0  (4-36) 

where  k  is  a  constant. 

Thus  the  slope  of  a  plot  of  log  (V/V0)  vs  log 
f(V/VQ)  immediately  yields  the  value  of  §,  and  moreover  de¬ 
parture  of  the  curve  from  linearity  indicates  the  range  of 
validity  of  the  fit.  A  sample  plot  is  shown  in  Fig. 4. 3  for 
aluminum,  the  constants  c  and  s  (Eq.4.35)  being  taken  from 
a  fit  to  the  data  of  Walsh  et  al.^  The  Mumaghan  expression 
Eq.4-33  is  a  good  fit  over  the  range  1  >  V/V0  >0.6,  i.e. 
over  a  much  wider  range  than  that  in  which  entropy  varia¬ 
tions  may  be  safely  neglected. 

Some  interesting  comparisons  are  possible. 

39 

Ryabinin  has  found  a  relationship  between  the  Murnaghan 
exponent  £  and  the  Gruneisen  coefficient  7  through  the 
Debye  theory 


$  -  2y  +  ^  (4.37) 

Using  7  (Table  II)  for  aluminum  leads  to  a 

value  of  £  *  4.33,  which  is  in  excellent  agreement  with  the 
value  £  ■*  4.34  from  Fig. 4. 3. 

40 

Borelius  has  also  connected  the  Gruneisen  ratio 
to  compressibility  data  by  thermodynamic  considerations, 
and  has  obtained  an  average  value  of  70  =  2,2  for  face- 
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centered  cubic  metals  and  -y0  ■  1.2  for  body- centered  cubic 
metals,  which  is  in  fair  agreement  with  values  of  70  ob¬ 
tained  by  Valsh  et  al. ^  One,  therefore,  expects  that  the 
Murnaghan  exponent  will  be  near  £  m  4.7  for  face-centered 
cubic  metals  and  £  -  2.7  for  body-centered  cubic  metals. 
This  leads  to  some  interesting  possibilities  in  scaling 
behavior  in  one  material  to  that  in  another  of  the  same 
crystal  structure.  An  application  will  be  mentioned  in 
Section  VII. 
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SECTION  V 


EXAMPLES  OF  WAVE- PROPAGATION  CALCULATIONS 


In  the  previous  three  sections,  a  discussion  was 
given  of  two  computational  methods,  viz.  the  finite  differ¬ 
ence  method  and  the  characteristic  method,  and^three  equa¬ 
tions  of  state,  viz.  the  Mie -Gruneisen ^equat ion,  the  e las tic - 


plastic  equation,  and  th6  Murnaghan  equation^  Solution  of  a 
particular  problem  requir&s  ehuice  ulr'the" ffiethod  of  computa¬ 
tion,  and  equation  of  state  best  suited  to  the  conditions  of 
the  problem.  In  this  section,  a  few  example  calculations 
will  be  described  in  some  detail,  which  will  illustrate  some 
of  the  advantages  and  disadvantages  of  each  method  of  solu¬ 
tion.  No  attempt  will  be  made  to  lay  down  general  rules  for 
choice  of  computational  method  or  equation  of  state,  however, 
and  each  problem  must  be  considered  separately  in  this  regard. 


Computations  were  performed  for  the  following  plate 


impact  problem: 

Driver  Plate  Thickness, 
Driver  Plate  Velocity, 
Driver  Plate  Material, 
Target  Plate  Thickness, 
Target  Plate  Material, 


0.9527  cm 

1.262  x  104  cm/ sec  *i  ti  FfS 

Aluminum 

1.2634  cm  //f? 

Aluminum 


The  expected  peak  stress  is  very  low,  about  10 
k i lobar s,  and  the  plate  thicknesses  are  such  that  rarefactions 
will  not  be  able  to  attenuate  the  shock  initially  propagating 
into  the  target  plate.  Stress  profiles  at  various  times  were 
required  for  several  reverberations  of  the  waves  in  the  tar¬ 
get  plate. 


Due  to  the  many  elastic-plastic  wave  interactions 
which  occur  on  reflection  of  the  wave  systems  at  free  sur- 
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faces,  and  on  crossing  of  wave  systems,  the  characteristic 
method  was  expected  to  be  very  complex,  and  not  suited  to 
automatic  computation.  A  semimanual  method  was  used,  as 
described  in  Section  IV.  The  elastic  and  plastic  portions 
of  the  equation  of  state  (Fig. 4.1)  could  be  approximated  by 
straight  line  segments,  due  to  the  very  small  curvature  of 
these  portions  over  the  small  compression  range  encompassed 
by  this  problem.  Nevertheless,  the  characteristic  mesh  con¬ 
sisted  of  more  than  200  wave  interactions.  Some  of  the  char¬ 
acteristics  are  shown  in  Fig. 5.1. 

Corresponding  stress-Lagrangian  position  profiles 
at  several  times  after  impact  are  shown  in  Fig. 5. 2.  Also 
shown  are  profiles  computed  by  the  finite  difference  method 
using  several  mesh  sizes.  Finally  for  comparison,  wave  pro¬ 
files  computed  by  the  finite  difference  method  using  the 
Murnaghan  equation,  i.e.,  material  with  zero  shear  strength, 
are  shown. 

Several  interesting  observations  may  be  made.  The 
elastic-plastic  interactions  lead  to  much  more  dispersion  of 
the  loading  and  unloading  waves  than  ir  a  material  of  zero 
shear  strength.  The  elastic-plastic  interactions  also  lead 
to  a  reduction  in  peak  stress  on  each  subsequent  wave  rever¬ 
beration.  Both  the  wave  dispersion  and  attenuation  are  re¬ 
produced  in  the  finite  difference  solutions.  The  fine  de¬ 
tails  of  the  waves  are  not  reproduced,  the  finite  difference 
solutions  tending  to  smooth  the  discontinuities.  Neverthe¬ 
less,  the  gross  features  of  the  wave  profiles  are  reproduced 
to  a  surprising  degree,  even  by  the  calculation  with  the 
coarsest  mesh  size.  Thus  the  peak  stress,  wavelength  and 
total  momentum  are  all  approximated  quite  closely.  The  de¬ 
tails  of  the  elastic-plastic  waves  are  reproduced  to  dif¬ 
ferent  degrees  by  calculations  with  different  mesh  sizes. 
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Only  the  calculation  with  the  finest  mesh  size  (350  meshes 
total)  was  able  to  approximate  the  main  elastic-plastic  dis¬ 
continuities  reasonably  well. 

In  choosing  the  calculation  method  and  mesh  size, 
it  is  necessary  to  decide  on  the  detail  desired  in  the  solu¬ 
tion.  The  characteristic  solution  gives  the  most  detailed 
representation  but  is  extremely  cumbersome  and  tedious  in 
the  present  case.  It  is  doubtful  if  the  fine  detail  shown 
by  the  characteristic  solution  would  actually  be  observed  in 
a  polycrystalline  material  where  a  certain  amount  of  disper¬ 
sion  and  smoothing  may  be  expected  due  to  the  presence  of 
grain  boundaries,  different  crystal  orientations,  natural 
viscosity,  etc.  For  the  finite  difference  method,  the  com¬ 
putation  time  increases  approximately  as  the  inverse  of  the 
square  of  the  mesh  size,  if  rezoning  or  the  time  saving  de¬ 
vices  are  not  used.  Thus,  for  a  mesh  size  Ar  =  0.0063  cm, 
the  computation  time,  using  an  I.B.M.  7090  digital  computer, 
was  about  40  minutes.  For  a  mesh  size  Ar  *  0.0126  cm,  the 
computation  time  was  about  11  minutes.  For  a  mesh  size 
Ar  =  0.063  cm,  the  computation  time  was  only  1.8  minutes, 
most  of  which  was  taken  up  in  input  and  output  routines. 

Neglecting  the  shear  strength  in  the  present  prob¬ 
lem  leads  to  rather  considerable  error,  since  both  the  peak 
stress  induced  by  the  impact  and  the  wave  velocities  are  in¬ 
correct.  At  later  times,  the  wave  profiles  bear  little  re¬ 
semblance  to  the  correct  solution.  Even  the  calculation 
with  the  coarsest  mesh  (35  meshes  total)  but  with  the 
elastic-plastic  equation  of  state,  gives  a  much  better  repre¬ 
sentation  of  the  wave  profiles. 

Another  computation  was  performed  for  the  following 
plate  impact  problem: 
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Driver  Plate  Thickness, 
Driver  Plate  Velocity, 
Driver  Plate  Material, 
Target  Plate  Thickness, 
Target  Plate  Material, 


0.3175  cm 

1.63  x  10*  cm/ sec 

Aluminum 

12  cm 

Aluminum 


In  this  case,  the  higher  driver  plate  velocity  is 
expected  to  lead  to  a  peak  stress  of  about  146  kilobars. 
Entropy  variations  are  expected  to  be  small.  On  the  other 
hand,  the  peak  stress  is  much  larger  than  the  yield  stress. 
Therefore  the  Murnaghan  equation  is  likely  to  be  adequate. 

In  this  particular  problem,  shock  attenuation  by  the  fol¬ 
lowing  rarefaction  was  of  interest,  and  successive  reflec¬ 
tion  at  the  target  plate  surfaces  was  not  required.  The 

characteristic  method  is  ideally  suited  to  this  problem, 

21 

and  Fowles  has  given  a  simple  analytical  solution  neglecting 
entropy  variations. 

Comparisons  of  Fowles1  solution  with  three  finite 
difference  solutions  are  shown  in  Fig. 5. 3  where  pressure- 
Eulerian  position  profiles  are  shown  for  six  times  after 
impact.  The  solutions  were  obtained  for  mesh  sizes  Ar  = 
0.0397;  0.0794;  and  0.1058  cm,  for  which  the  driver  plate 
contained  8,  4,  and  3  meshes,  respectively. 

Some  of  the  features  inherent  in  the  finite  dif¬ 
ference  method  are  illustrated  by  Fig. 5. 3.  The  shock  is 
spread  over  about  four  meshes,  and  therefore  the  head  of  the 
shock  is  retarded,  compared  to  the  exact  solution.  The  rare¬ 
faction  is  also  distorted,  the  tail  being  retarded,  although 
the  artificial  viscosity  was  set  zero  on  rarefactions*  The 
net  effect  is  that  the  peak  pressure  is  attenuated  more 
rapidly  than  in  the  correct  solution.  Partially  off- setting 
this  effect  is  the  fact  that  the  finite  difference  solution 
overshoots  the  peak  pressure  at  the  shock,  and  introduces  a 
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series  of  oscillations  into  the  subsequent  pressure  profile. 

Inherent  in  the  finite  difference  solution,  there¬ 
fore,  is  an  uncertainty  in  the  peak  pressure  and  pulse  width, 
the  uncertainty  increasing  with  mesh  size.  Nevertheless  the 
gross  wave  shape  is  reproduced  quite  well.  Calculations  of 
momentum  distribution  are  therefore  likely  to  be  handled 
quite  well  by  the  finite  difference  method,  but  calculations 
involving  fracture  criteria  are  much  more  sensitive  to  the 
peak  pressure  and  pulse  width,  and  the  finite  difference 
method  may  not  be  adequate  if  accuracy  is  desired. 

It  should  be  pointed  out  that  much  greater  accuracy 
is  necessary  when  the  peak  stress  is  attenuated  by  a  large 
amount.  For  example,  if  the  initial  induced  stress  at  the 
interface  is  100  kilobars,  and  the  uncertainties  introduced 
due  to  inaccuracies  in  the  equation  of  state  and  by  the 
finite  difference  method  are  five  per  cent,  then  if  the  peak 
stress  is  subsequently  attenuated  to  10  kilobars,  the  uncer¬ 
tainty  in  the  result  will  be  of  the  order  of  50  per  cent. 

In  calculations  of  fracture  thresholds  under  such  conditions, 
therefore,  it  is  necessary  to  ensure  great  accuracy  if  the 
results  are  to  be  meaningful,  and  characteristic  methods  are 
preferred. 

Introduction  of  shear  strength  into  the  equation 
of  state  does  not  lead  to  any  significant  changes  in  the 
finite  difference  solutions  in  this  case.  Only  small  dif¬ 
ferences  would  be  introduced  into  the  characteristic  solu¬ 
tion.  The  finite  difference  solutions  with  the  above  mesh 
sizes  are  not  capable  of  reproducing  these  differences. 

Other  examples  of  wave  propagation  calculations 
are  included  in  the  next  two  sections. 
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SECTION  VI 


SPALL  FRACTURE 


One  of  the  most  important  aspects  of  the  response 
of  solid  bodies  to  intense  impulsive  loads  is  that  internal 
fractures  may  be  produced.  In  some  cases  where  the  frac¬ 
tures  are  close  to  a  surface,  pieces  of  the  solid  may  be 
tom  loose  and  fly  off  at  considerable  speed.  This  phe¬ 
nomenon  is  broadly  termed  spallation. 

In  order  to  predict  the  occurrence  of  spall  frac¬ 
ture,  it  is  necessary  to  know  the  time  dependent  stress  dis¬ 
tribution  through  the  solid  body.  This  may  be  computed 
from  the  loading  and  boundary  conditions  once  an  adequate 
constitutive  equation  has  been  found  for  the  particular 
material  under  study.  It  is  also  necessary  to  know  the  frac¬ 
ture  criterion  for  the  material,  i.e.,  conditions  under  which 
mechanical  failure  occurs. 

This  section  will  be  devoted  to  a  brief  discussion 
of  fracture  criteria,  and  their  evaluation. 

6.1  Fracture  Criteria 

The  most  commonly  used  fracture  criterion  considers 
that  failure  occurs  when  the  greatest  normal  tensile  stress 
reaches  a  critical  value.  This  is  by  no  means  the  only 
fracture  criterion  in  use.  Other  criteria  are  based  on  criti¬ 
cal  values  of  shear  stress,  relative  elongation,  or  work  done 
in  deformation.  If  attention  is  restricted  to  one  particular 
configuration,  in  this  case  uniaxial  strain,  then  any  one  of 
these  criteria  may  be  chosen  providing  that  the  criterion  is 
evaluated  under  identical  conditions  of  uniaxial  strain.  The 
most  convenient  parameter  in  this  case  is  the  normal  stress. 
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While  this  leads  to  a  self-consistent  fracture  criterion 
under  conditions  of  uniaxial  strain,  it  cannot  necessarily 
be  related  directly  to  fracture  criteria  obtained  in  other 
configurations,  e.g.,  uniaxial  stress. 

The  most  convenient  experiment  for  evaluating  the 
fracture  criterion  under  conditions  of  uniaxial  strain  is 
the  plate  impact  experiment.  A  flat  target  plate  of  the 
material  under  study,  of  lateral  extent  large  compared  to  its 
thickness,  is  impacted  face  on  by  a  thin,  flat  driver  plate. 
Uniaxial  strain  is  maintained  near  the  center  of  the  plates 
for  a  sufficient  time  to  produce  tension  stresses  in  the  in¬ 
terior  of  the  target  plate.  The  magnitude  and  duration  of 
the  tension  stresses  may  be  controlled  by  choice  of  driver 
plate  velocity  and  plate  thicknesses. 

A  number  of  such  experiments  have  been  performed 

at  various  laboratories.  Figure  6.1  shows  sections  through 

1/4-inch  thick  copper  target  plates  impacted  at  increasing 

velocities  by  1/8- inch  thick  copper  driver  plates  by 

22 

Lundergan  et  al.  Great  care  was  exercised  in  these  tests 
to  ensure  that  the  impact  was  plane  and  that  air  between  the 
plates  was  removed  (See  Section  VII). 

At  a  driver  plate  velocity  of  0.52  x  104  cny'sec, 
only  a  few  isolated  cracks  appear  in  the  expected  spall 
plane.  At  a  driver  plate  velocity  of  0.896  x  104  cm/ sec 
the  cracking  is  much  more  extensive.  It  is  interesting  to 
observe  that  not  all  of  the  cracks  lie  on  a  sharply  defined 
plane,  and  do  not  all  join  to  form  a  smooth  fracture  surface. 
Rather,  cracks  lie  on  different  planes  and  tend  to  overlap, 
leaving  small  tongues  of  material  connecting  the  spall  to 
the  remaining  target.  At  higher  driver  plate  velocities, 
the  cracking  is  apparently  more  severe,  until  at  a  driver 
plate  velocity  above  1.141  x  104  cny^sec  the  spall  is  completely 


ASD-TDR-62-399 


50 


severed  from  the  remaining  target. 

From  the  photomicrographs  (Fig.6.1),  it  is  clear 
that  the  target  sustains  increasing  damage  as  the  driver 
plate  velocity,  and  hence  the  induced  stress,  is  increased. 

It  is  very  difficult  to  assess  the  degree  of  damage  from 
the  photomicrographs.  It  was  felt  that  the  residual  strength 
of  the  target  on  the  spall  plane  would  give  a  good  indication 
of  the  degree  of  damage.  Consequently,  small  tension  speci¬ 
mens  were  cut  from  spalled  target  plates  kindly  supplied  by 
Dr .C.D.Lundergan  of  Sandia  Corporation,  and  the  stress  neces¬ 
sary  to  separate  the  spall  from  the  target  was  measured  in 
each  case.  Specimens  were  mounted  in  special  jaws  in  a 
tensile  testing  machine  as  shown  in  Fig. 6. 2. 

Results  are  shown  in  Fig. 6. 3.  Data  from  specimens 
cut  from  the  same  target  showed  considerable  scatter  with  a 
tendency  for  the  residual  strength  to  be  somewhat  higher 
near  the  edges  of  the  spalled  region  than  in  the  center. 

Some  uncertainty  may  have  resulted  from  the  fact  that  ferric 
chloride  etchant  had  been  introduced  into  the  cracks  during 
the  metallurgical  examination,  and  the  etchant  may  have  re¬ 
duced  the  residual  strength  of  the  specimens.  However,  it 
was  felt  that  the  effect  of  the  etchant  was  small,  and  more 
or  less  uniform  in  all  the  specimens,  and  did  not  invalidate 
the  general  trend  of  the  data. 

The  trend  of  the  data  is  clear.  For  driver  plate 
velocities  below  about  0.5  x  10*  cny'sec  no  damage  is  evident 
on  impacted  target  plates,  and  the  strength  should  be  that  of 
the  virgin  material.  Increasing  damage  is  sustained  with  in¬ 
creasing  driver  plate  velocity  up  to  a  velocity  of  about 
0.8  x  10*  cny'sec.  Beyond  this  velocity,  no  further  change 
is  apparent  in  the  residual  strength,  and  the  only  observable 
effect  is  a  slight  widening  of  the  cracks  (Fig. 6.1),  up  to  a 
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driver  plate  velocity  of  1.41  x  104  cn/sec,  at  which  velocity 
the  spall  is  completely  detached.  Corresponding  thresholds 
may  be  arbitrarily  defined:  these  will  be  termed  the  lower 
fracture  threshold,  below  which  no  damage  is  apparent  in 
the  target;  the  upper  fracture  threshold,  above  which  the 
cracks  are  completely  formed;  and  the  spall  separation 
threshold,  above  which  the  spall  separates.  For  driver 
plate  velocities  above  the  spall  detachment  threshold,  the 
velocity  of  the  ejected  spall  is  roughly  linearly  related  to 
the  driver  plate  velocity  (Fig. 6. 4). 

In  order  to  study  the  fracture  phenomenon  further, 
it  is  necessary  to  inquire  into  the  stress- time  history  at 
the  fracture  location.  Unfortunately,  the  Hugoniot  (Section 
IV)  for  this  material  has  not  been  measured  in  the  low  stress 
region.  It  is  only  possible  to  make  rough  estimates  of  the 
stresses  with  the  aid  of  a  constitutive  equation  constructed 
from  the  dilatational  adiabat  from  Section  4.2,  and  a  hand¬ 
book  value  of  the  yield  stress.  While  this  will  not  lead  to 
very  accurate  quantitative  data  on  the  fracture  stress,  it 
can  serve  to  illustrate  qualitative  aspects  of  the  phenomenon. 

The  characteristic  net  for  a  typical  experiment  will 
appear  as  in  Fig. 6. 5,  where  it  is  assumed  that  no  fracture 
has  occurred.  (Some  of  the  characteristics  have  been  omitted 
for  clarity.)  Corresponding  plots  of  the  stress-time  history 
at  locations  at  the  center  of  the  target  plate,  and  one 
quarter  of  the  target  plate  thickness  from  the  interface  have 
been  plotted  in  Fig. 6. 6.  Qualitatively,  the  stress  time  pro¬ 
files  are  nearly  square  waves.  The  magnitude  of  the  tensile 
stress  is  directly  related  to  the  driver  plate  velocity. 

The  duration  of  the  first  tensile  cycle  is  directly  pro¬ 
portional  to  the  distance  of  the  location  from  the  interface 
(or  free  back  surface),  provided  that  no  fractures  occur. 
Moreover,  as  long  as  no  fractures  occur,  the  positive  and 
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negative  cycles  will  repeat  with  only  minor  modification 
until  waves  from  the  edges  of  the  plates  destroy  the  uni¬ 
axial  strain  configuration  and/or  the  target  plate  strikes 
some  object  in  its  flight  path. 

Using  the  approximate  equation  of  state  described 
above,  the  magnitude  of  the  stress  in  the  first  tension 
cycle  corresponding  to  the  lower  fracture  threshold  is  about 
8.5  kilobars,  and  corresponding  to  the  upper  fracture  thres¬ 
hold  is  about  14  kilobars. 

Despite  the  apparent  simplicity  of  the  test  con¬ 
figuration,  interpretation  of  the  observed  phenomena  is  not 
at  all  straightforward.  The  fact  that  the  target  sustains 
increasing  damage  in  the  range  between  lower  and  upper  frac¬ 
ture  thresholds,  as  defined  above,  suggests  that  the  fracture 
phenomenon  is  time  dependent.  Such  time  dependence  is,  of 
course,  familiar  from  tensile  tests  carried  out  by  many 
experimenters,  'over  time  ranges  varying  from  milliseconds 
to  months.  Simple  reaction  rate  kinetics  has  been  applied  to 
such  relatively  long  time  fracture  data,  as  well  as  creep 
data,  with  considerable  success. 

By  considering  the  probability  of  bond  rupture 
under  an  applied  load,  it  is  possible  to  deduce  an  expression 

for  the  fracture  time  t  ,  following  arguments  such  as  those  of 

27  c 

Beuche,  of  general  form 


where  to  is  the  vibrational  frequency  of  the  crystal  lattice, 
and  is  thus  the  number  of  times  per  second  that  a  particular 
bond  reaches  a  condition  under  which  rupture  may  occur,  k  is 
the  Bo  Italian  constant,  T  is  the  temperature,  and  t  is  the 
magnitude  of  the  activation  barrier,  which  is  a  function  of 
the  applied  stress  o. 
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Many  experimenters,  including  Zhurkov  et  al, 
found  that  experimental  creep  and  fracture  data  could  be 
fitted  by  an  expression 

tc  -  i  «p  <6-2> 

28 

Other  experimenters,  including  Dorn,  noted  that  creep  data 
could  be  fitted  by  an  expression  equivalent  to 

tc  “Se*P  -p’°>  (6-3) 

i.e.,  with  the  constant  0  of  Eq.6.2  proportional  to  tempera¬ 
ture.  Measured  values  of  range  from  the  activation  en- 

29 

ergy  of  sublimation  to  that  of  self-diffusion.  Anderson 
in  a  discussion  of  some  of  the  apparent  discrepancies  in  the 
literature  has  postulated  a  dual  mechanism,  self -diffusion 
dominating  at  low  rates  and  high  temperatures,  when  Eq.6.3 
may  be  applicable,  and  sublimation  dominating  at  high  rates, 
with  Eq.6.2  as  the  governing  equation.  Without  entering 
Into  a  discussion  of  these  points,  it  might  be  noted  that 
under  brittle  fracture  conditions  found  in  spallation,  the 
sublimation  process  is  likely  to  be  the  dominant  one. 
Moreover,  all  of  the  Sandia  experiments  quoted  above  in¬ 
volve  low  stresses,  so  that  the  local  temperature  during 
fracture  is  very  slightly  below  ambient  temperature  in 
all  cases,  in  which  case  Eq.6.2  and  6.3  are  equivalent. 

Extrapolation  of  Eq.6.2  below  0(10  3)  seconds  is 
not  necessarily  valid.  Observed  cracks  under  spall  condi¬ 
tions  (Fig. 6.1)  are  transgranular  and  appear  to  run  continu¬ 
ously  for  10  to  20  grain  diameters.  For  very  short  fracture 
times,  the  crack  propagation  speeds  must  be  very  high,  near 
the  sonic  velocity  of  the  material.  Waldorf has  suggested 
that  the  stress  necessary  to  propagate  cracks  at  such  ve- 
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locities  must  be  higher  than  that  given  by  Eq.6.2,  although 
the  rate  of  increase  proposed  by  Waldorf  seems  too  high. 
Valluri  and  Williams,  1  in  a  discussion  of  spall  fracture, 
point  out  that  the  lower  limit  on  the  fracture  time  must  be 
of  the  order  of  the  vibrational  period  of  the  crystal  lat¬ 
tice,  i.e.  0(10-13)  seconds,  and  the  corresponding  stress 
must  be  the  theoretical  strength  of  the  material,  in  this 
case  about  200  kilobars.  Experimental  evidence  in  support  of 
such  a  limit  has  recently  been  obtained  by  McQueen  and 
Marsh, ^  who  have  obtained  a  limiting  tensile  strength  of 
233  kilobars.  This  stress  is  again  much  higher  than  that  pre¬ 
dicted  by  Eq.6.2. 

Unfortunately,  relatively  long  time  creep  and 
rupture  data  is  not  available  for  the  material  under  con¬ 
sideration,  i.e.  half-hard  E.T.P.  copper.  However,  the 
straight  line  in  Fig. 6. 7  was  plotted  from  Eq.6.2,  using  the 
sublimation  energy,  the  lattice  vibrational  frequency,  and 
a  value  of  0  such  that  the  curve  fitted  the  "static"  ultimate 
stress  (handbook  value)  in  the  range  1  to  102  seconds.  The 
resultant  value  of  P  is  of  correct  order  as  predicted  by 
consideration  of  ionic  volumes.  The  above  considerations 
then  lead  to  the  expectation  that  the  qualitative  behavior 
at  shorter  times  will  be  roughly  as  shown  by  the  dotted  ex¬ 
tension  in  Fig. 6. 7. 

Clearly  Fig. 6. 7  is  merely  qualitatively  correct, 
and  is  intended  only  to  further  the  discussion  of  the  frac¬ 
ture  mechanism.  The  behavior  suggested  by  Fig. 6. 7  would 
lead  to  a  reasonably  consistent  explanation  of  the  phenomena 
observed  in  spallation  experiments. 

For  stresses  below  the  lower  fracture  threshold, 
i.e.  below  about  8  kilobars,  the  time  required  for  fracture 
would  be  in  excess  of  0(10  3)  seconds.  In  the  experiment, 
rarefaction  waves  from  the  edges  of  the  plates  attenuate 
stresses  in  the  target  in  times  of  order  10  5  seconds.  Thus 

no  damage  would  be  noticeable  in  the  target. 
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On  the  other  hand,  if  the  tensile  stress  is  in  ex¬ 
cess  of  the  upper  fracture  threshold,  i.e.  about  14  ki lobar s, 
the  fracture  time  would  be  less  than  0(10  8)  seconds,  which 
is  very  short  compared  to  the  duration  of  the  first  tensile 
cycle,  i.e.  0(10“®)  seconds.  Thus  to  all  intents  and  pur¬ 
poses,  fracture  would  occur  instantaneously  on  the  time 
scale  of  the  experiment,  and  the  fracture  time  can  be  ignored. 

In  the  region  between  the  lower  and  upper  fracture 
thresholds,  however,  the  fracture  time  would  be  significant, 
and  would  lead  to  just  the  type  of  behavior  observed,  i.e. 
increasing  damage  for  increasing  stress  level.  Interpreta¬ 
tion  of  the  fracture  data  in  this  region  is  very  difficult 
because  the  loading  is  now  no  longer  simple.  As  the  cracks 
grow  they  will  originate  release  waves  which  affect  the  stress 
wave  propagation  in  the  target,  and  the  wave  profiles  shown  in 
Fig. 6. 3  will  be  profoundly  altered.  There  is  also  the  pos¬ 
sibility  that  subsequent  reverberations  of  the  stress  waves 
would  lead  to  further  growth  of  the  fractures. 

The  effect  of  loading  time  on  the  fracture  forma¬ 
tion  process  should  be  observable  in  the  region  between  lower 
and  upper  fracture  threshold,  and  various  experiments  have 
been  proposed  to  observe  the  effect.  The  duration  of  the 
first  tensile  cycle  can  be  varied  by  changing  the  geometrical 
scale  of  the  experiment.  Due  to  experimental  difficulties, 
the  scale  could  be  conveniently  altered  at  most  by  a  factor 
of  200,  i.e.  first  tensile  cycle  durations  from  about  2  x  10*e 
seconds  to  5  x  10~®  seconds.  If  Fig. 6. 7  is  reasonably  cor¬ 
rect,  the  upper  fracture  threshold  would  alter  by  only  10 
or  20  per  cent  in  this  range.  Thus,  determination  of  upper 
fracture  threshold  by  measurements  of  residual  strength 
(Fig. 6. 3)  would  have  to  be  carried  out  very  carefully,  if  a 
fracture  time  effect  is  to  be  observed.  Conversely,  if  a 
fracture  time  effect  is  not  observed  in  relatively  crude  ex- 
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periments,  this  would  serve  as  verification  that  Fig. 6. 7  is 
qualitatively  correct,  i.e.  that  the  slope  in  the  vicinity 

10"8  to  10~°  seconds  is  not  very  much  greater  than  shown. 

22  32 

Preliminary  experiments  *  have  failed  to  detect  a  time  de¬ 
pendence. 

There  is  an  important  difference  between  the  rela¬ 
tively  long  time  tensile  tests,  results  of  which  lead  to  ex¬ 
pressions  of  the  form  of  Eq.6.2,  and  the  spall  experiments. 
The  former  experiments  are  usually  carried  out  tinder  dead 
weight  loading.  Thus,  once  a  mechanical  instability  has 

formed,  the  test  is  carried  out  to  complete  failure  under 

25 

constant  load,  or,  as  in  the  tests  of  Zhurkov,  under  a 
load  programmed  to  yield  a  constant  nominal  true  stress,  de¬ 
fined  as  the  load  divided  by  the  neck  area. 

On  the  other  hand,  in  the  spall  experiment,  as  the 
fractures  develop,  widening  of  the  cracks  may  be  expected  to 
initiate  release  waves  which  will  tend  to  attenuate  the 
nominal  stress.  One  thus  expects  a  difference  in  behavior 
during  the  last  stages  of  the  fracture  process.  Above  the 
upper  fracture  threshold,  the  cracks  may  be  expected  to 
grow  very  rapidly,  but  widening  of  the  cracks  may  be  expected 
to  relieve  the  nominal  stress  almost  immediately.  Since 
cracks  are  initiated  at  suitable  nuclei  dispersed  randomly 
in  the  material,  and  their  propagation  directions  are  influ¬ 
enced  by  local  grain  structure,  the  cracks  are  not  all  co- 
planar,  and  do  not  all  join  smoothly,  and  the  small  inter¬ 
connecting  tongues  observed  in  Fig. 6.1  remain.  At  this 
point  the  spall  is  moving  at  nearly  driver  plate  velocity, 
while  the  remaining  target  is  at  nearly  zero  velocity. 
Metallurgical  observations  reveal  that  the  tongues  appear 
to  have  been  drawn  out  plastically,  and  occasionally  evi¬ 
dences  of  local  melting  are  observed  at  the  tips  of  ruptured 
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tongues.  Thus  the  final  stage  of  the  fracture  process  ap¬ 
pears  to  involve  a  "quasi- static"  plastic  extension  of  the 
tongues ,  and,  when  the  imparted  impulse  is  sufficiently 
high,  a  ductile  rupture  of  the  tongues  occurs. 

Complete  separation  of  the  spall  from  the  target 
is  generally  complicated  by  edge  effects.  Figure  6.8  shows 
a  section  near  the  edge  of  the  target  plate.  In  the  Sandia 
series  of  tests,  a  groove  was  cut  to  the  expected  depth  of 
the  spall  in  order  to  delimit  the  test  area.  Since  the 
groove  introduces  a  lateral  free  surface,  which  must  remain 
at  zero  stress,  the  maximum  tension  stress,  and  hence  the 
fractures,  are  not  developed  adjacent  to  the  groove.  This 
ring  of  undamaged  material  must  thus  be  sheared  if  the  spall 
is  to  be  detached. 

When  no  groove  is  provided,  the  cracks  are  limited 
to  the  area  of  the  target  plate  which  is  subjected  to  the 
tension  stress,  and  in  order  to  detach  the  spall  it  is  again 
necessary  to  shear  material  between  the  edge  of  the  cracked 
region  and  the  surface.  In  similar  tests  with  and  without 
grooves,  roughly  the  same  driver  plate  velocity  was  re¬ 
quired  to  detach  the  spall.  Tests  have  also  been  performed 
in  which  the  test  area  was  delimited  by  a  circumferential 
cut  clear  through  the  target  plate.  The  test  specimen  was  a 
good  fit  in  the  surrounding  plate.  In  this  case,  cracks  are 
formed  right  up  to  the  edge  of  the  test  piece  and  edge  ef¬ 
fects  are  eliminated.  A  reduction  in  spall  detachment  thres- 

22 

hold  of  15  per  cent  was  observed. 

It  should  be  noted  that  the  effect  of  the  edge 
shear  ring  is  relatively  more  important  when  the  spall  is 
thick,  and  when  the  spalled  area  is  small.  One  would  thus 
expect  a  size  effect  on  spall  detachment  threshold.  Such  a 
size  effect  has  in  fact  been  observed. 
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While  the  lower  and  upper  fracture  thresholds  may 
be  defined  in  terms  of  the  maximum  developed  tensile  stress, 
and  appear  to  be  only  weakly  dependent  on  stress  duration, 
the  spall  separation  threshold  will  be  a  function  of  the  target 
geometry  and  loading  configuration.  Not  only  will  the  spall 
separation  threshold  be  a  function  of  spall  size,  and  edge 
restraint,  but  also  of  the  number  and  strength  of  intercon¬ 
necting  tongues.  Thus  the  grain  structure,  which  determines 
the  crack  configuration  should  be  very  important.  However, 
equally  important  will  be  the  planarity  of  the  stress  wave 
system  in  the  target  plate.  If  the  waves  in  the  target 
plate  are  very  perturbed,  for  example  because  of  a  rough  or 
buckled  driver  plate  or  a  roughened  rear  surface  of  the 
target,  cracks  would  be  nucleated  at  more  widely  separated 
points,  and  a  much  higher  impulse  would  be  required  to 
separate  the  spall. 

For  stresses  beyond  the  upper  fracture  threshold, 
it  is  unlikely  that  the  tension  stress,  calculated  by  as¬ 
suming  that  no  fractures  occur,  would  actually  be  experienced 
by  the  material  on  the  spall  plane.  It  is,  therefore,  not 
correct  to  speak  of  a  spall  separation  threshold  stress. 

Rather,  based  on  the  above  hypothesis,  the  separation 
mechanism  depends  on  the  relative  average  velocities  of  the 
spall  and  remaining  target.  These  in  turn  depend  on  the 
stress  wave  profiles  in  the  spall  and  target  at  the  time  the 
fractures  occur.  A  single  quantity  which  is  adequate  to 
serve  as  a  measure  of  the  spall  separation  threshold  has  so 
far  not  been  defined,  but  it  is  possible  that  the  relative 
average  velocity  of  spall  and  target  may  prove  suitable  for 
this  purpose. 

The  spall  velocity,  when  separation  occurs,  may 
also  be  expected  to  be  a  function  of  the  edge  restraint  and 
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number  and  strength  of  interconnecting  tongues.  In  experi¬ 
ments  in  which  no  delimiting  grooves  or  rings  are  provided 
there  is  considerable  scatter  in  the  measured  spall  velocity, 
presumably  due  to  random  differences  in  edge  restraint.  In 
the  Sandia  series  of  experiments,  the  scatter  is  not  too 
great  (Fig. 6. 4).  Extrapolation  of  the  curve  to  the  abscissa 
indicates  that  considerable  energy  is  expended  in  forming 
the  completed  fracture.  However,  insufficient  data  is  avail¬ 
able  to  determine  the  amount  of  energy  expended  in  the  vari¬ 
ous  modes,  i.e.  formation  of  the  cracks,  rupture  of  the 
tongues,  shearing  of  the  edge  restraint,  and  contribution 
to  the  remaining  target  kinetic  energy.  More  work  is  neces¬ 
sary  to  determine  the  energy  balance  before  rational  pre¬ 
dictions  of  spall  velocity  can  be  made. 

6.2  Summary 

The  preceding  discussion  has  necessarily  been 
qualitative  and  speculative,  due  to  the  lack  of  careful  ex¬ 
perimental  data.  Nevertheless  a  hypothesis  has  been  ad¬ 
vanced  which  is  in  accord  with  most  of  the  experimental  ob¬ 
servations.  The  most  important  deductions  of  engineering 
importance  are  the  existence  of  three  thresholds;  i.e.,  the 
lower  and  upper  critical  fracture  threshold,  and  the  spall 
separation  threshold.  Corresponding  stresses  may  be  as¬ 
signed  to  the  first  two  thresholds,  but  not  to  the  last. 

1.  For  stresses  below  the  lower  critical  fracture 
threshold,  no  effect  on  the  target  material  is  evident,  either 
under  low  powered  metallurgical  examination,  or  in  load 
carrying  capacity.  This  threshold  is  the  simplest  to  de¬ 
termine  experimentally. 

2.  For  conditions  above  the  upper  critical  frac¬ 
ture  threshold,  fully  formed  cracks  are  observed  on  the 
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spall  plane,  but  the  spall  is  not  necessarily  detached  from 
the  target.  The  residual  load  carrying  capacity  of  the 
spall  plane  depends  on  the  grain  structure  of  the  material 
and  the  planarity  of  the  stress  waves  which  gave  rise  to  the 
cracks,  and  may  vary  very  widely.  The  upper  fracture  thres¬ 
hold  should  be  weakly  dependent  on  loading  time,  and  thus  on 
the  geometrical  scale  ot  the  experiment,  but  tor  engineering 
purposes,  it  is  likely  that  this  variation  can  be  neglected. 

3.  The  spall  separation  threshold  defines  condi¬ 
tions  under  which  the  spall  separates  from  the  target.  This 
threshold  is  of  prime  engineering  importance,  but  is  the 
most  difficult  to  determine.  A  single  quantity  to  serve  as 
a  measure  of  the  spall  separation  threshold  has  so  far  not 
been  found.  The  spall  separation  threshold  will  be  a  func¬ 
tion  of: 

a)  the  relative  average  velocity  of  spall  and  re¬ 
maining  target  at  the  time  of  crack  formation,  and  thus  of 
the  stress  wave  profiles. 

b)  the  strength  of  interconnecting  tongues  be¬ 
tween  adjacent  cracks  formed  by  the  stress  waves,  and  hence 
of  the  grain  structure,  and  planarity  of  the  stress  waves 
which  gave  rise  to  the  crack. 

c)  the  edge  restraints. 

Variations  of  more  than  an  order  of  magnitude  have  been  ob¬ 
served  in  this  threshold. 

4.  The  spall  velocity  will  be  a  function  of  the 
kinetic  energy  entrapped  in  the  spall,  and  hence  on  the 
wave  profile  at  the  time  that  the  fracture  is  completed. 

The  wave  profile  will  be  affected  by  the  fact  that  energy 
is  necessary  to  form  the  fracture. 
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SECTION  VII 


EFFECT  OF  AIR  BETWEEN  IMPACTING  FLAT  PLATES 

An  experimental  technique,  which  has  been  widely 
used  for  studying  wave  propagation  and  spall  fracture,  in¬ 
volves  impacting  a  flat  plate  of  the  material  under  study 
with  a  thin,  flat  driver  plate.  Many  such  experiments  have 
been  carried  out  in  which  the  driver  plate  is  propelled  by 
an  explosive  charge  or  exploding  foil,  without  provision 
for  removing  the  atmospheric  air  between  driver  plate  and 
target.  Since  the  air  is  progressively  compressed  as  the 
driver  plate  approaches  the  target,  it  might  be  expected 
that  the  wave  profile  induced  in  the  target  plate  will  dif¬ 
fer  significantly  from  that  induced  in  the  absence  of  air. 

It  is  desirable  to  investigate  this  effect. 

Evidently,  the  motion  of  the  air  will  be  very 
complex,  since  the  plates  are  of  finite  size,  and  the  gap 
around  the  edges  is  usually  unsealed.  The  resultant  three- 
dimensional  flow  would  be  extremely  difficult  to  analyze. 

In  order  to  fix  an  upper  limit  on  the  effect,  a 
simple  one -dimensional  analysis  was  undertaken.  This  would 
be  directly  applicable  if  the  edge  gap  were  sealed.  For  an 
unsealed  gap,  the  unsteady  air  motion  is  such  that  more  air 
will  flow  from  the  gap  near  the  edges  of  the  plate,  and  the 
one- dimensional  analysis  would  be  approached  only  near  the 
center  of  the  plate.  Near  the  edges,  conditions  would  more 
nearly  approach  conditions  with  no  air  present,  and  the  in¬ 
duced  wave  profiles  would  therefore  vary  from  the  center  to 
the  edges  of  the  plate,  leading  to  complex  three-dimensional 
stress  wave  profiles  in  the  target. 

Attention  was  focused  on  impact  in  aluminum  driver 
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and  target  plates,  with  expected  peak  pressures  of  the  order 
of  100  kilobars.  Thus  a  simple  equation  of  state  P  ■  f(V), 
fitted  to  the  experimental  Hugoniot  (Section  4.2)  was  used, 
and  wave  propagation  in  the  plates  was  analyzed  using  the 
method  of  characteristics  (Section  III).  While  the  analysis 
can  be  easily  extended  to  lower  stresses  by  using  an 
elastic-plastic  equation  of  state,  much  higher  pressures 
could  be  analyzed  only  if  entropy  changes  are  allowed, 
either  by  using  the  finite-difference  method,  or  by  extending 
the  characteristic  method  to  accommodate  entropy  variations. 

7.1  Perfect  Gas  Analysis 

In  order  to  investigate  the  probable  magnitude  of 
the  effect,  it  was  first  assumed  that  the  air  behaved  as  a 
perfect  gas  throughout  the  motion.  It  was  also  assumed 
that  the  driver  plate  started  into  motion  impulsively.  When 
explosive  propulsion  techniques  are  used,  the  driver  plate 
acceleration  is  very  high,  so  that  the  compression  wave  in 
the  air  ahead  of  the  driver  steepens  into  a  shock  very 
rapidly,  and  actual  conditions  probably  approach  the  idealiza¬ 
tion  after  only  a  few  millimeters  of  driver  plate  motion. 

At  the  velocities  investigated,  the  pressure  in  the 
air  behind  the  first  shock  is  quite  low.  It  was  assumed  that 
the  driver  plate  was  in  a  stress-free  condition  at  this  time; 
i.e.,  that  stress  waves  due  to  propulsion  have  damped  out 
before  the  pressure  in  the  entrapped  air  becomes  appreciable. 

The  shock  wave  preceding  the  driver  plate  will  re¬ 
flect  successively  at  the  target  and  driver  plates,  raising 

the  pressure  in  successive  steps,  and  inducing  stress  waves 

33 

in  the  plates  (Fig. 7.1).  Evans  and  Evans  have  shown  that 
under  conditions  of  successive  shock  reflection  between  ap¬ 
proaching  rigid  walls,  the  bulk  of  the  entropy  rise  occurs 
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in  the  first  four  shock  compressions,  and  that  the  subsequent 
motion  is  very  nearly  isentropic.  Since  the  pressure  in  the 
air  is  still  relatively  low  at  this  point,  the  problem  was 
simplified  as  follows.  The  first  four  shock  reflections 
were  treated  by  the  method  of  Evans  and  Evans,  assuming  the 
plates  to  be  rigid,  with  the  driver  plate  moving  at  con¬ 
stant  velocity.  In  the  subsequent  motion,  the  air  was  as¬ 
sumed  to  behave  isentropically,  and  the  transit  times  of 
waves  across  the  air  gap  were  neglected,  but  the  plates  were 

considered  compressible.  A  somewhat  similar  approximation 

3  ^ 

has  been  made  by  Winter,  and  has  been  found  to  give  good 
results  in  an  analysis  of  the  motion  of  a  rigid  piston  in 
an  air-filled  tube. 

Evans  and  Evans  give  the  following  relations  for 
the  pressure  P  and  specific  volume  V  behind  the  nth  shock 
relative  to  initial  conditions: 


where 


A 


it  ±  l 

TTx  -  1 


and  p. 


^  +  1 

7-l 


Here  P  is  the  density  of  the  air,  iri  is  the  pres¬ 
sure  ratio  across  the  first  shock  and  7  is  the  ratio  of  the 
specific  heats  of  the  air,  taken  to  be  1.4.  The  value  of 
tti  is  obtained  from  the  initial  driver  plate  velocity  by 
the  Rankine-Hugoniot  relations  (for  the  air  shock) .  Repre¬ 
sentative  results  are  listed  in  Table  VI  for  several 
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driver  plate  velocities* 

During  the  subsequent  motion,  the  air  is  assumed 
to  be  described  by  the  isentropic  relation 


_L.  -  ^  =  constant  (7*1) 

Py  P? 

where  (  )4  pertains  to  conditions  after  the  fourth  shock  re¬ 
flection. 

In  order  to  use  the  method  of  characteristics,  it 
is  now  only  necessary  to  modify  the  interface  condition  Eqs. 
3.12a  and  3.12b,  Section  III,  by  specifying  that  the  inter¬ 
face  pressure  is  given  by  Eq.7.1  and 


Here  R  is  the  Euler ian  position  of  the  surface, 

(  ) '  refers  to  the  right-hand  plate  and  2X0  is  the  initial 
plate  spacing  before  commencement  of  the  motion,  at  which 
time  the  initial  (ambient)  air  density  is  P0* 

The  resultant  equations,  one  each  of  Eqs. 3. 8,  3.9, 
and  3.13  for  each  plate,  together  with  Eqs. 3. 10,  3.11,  7.1, 
and  7.2  are  implicit.  The  simplest  practical  method  of  solu¬ 
tion  is  to  estimate  the  positions  of  the  surfaces  by  modifying 
Eqs. 3. 9  and  3.13  to 


r  —  r 

num  m-i jm  ^ 

t  _  ”  t  ,  m-i,m 

in,®  ra-i,ra 


(7.3) 


R-R  +  U  (C-t  )  (7.4) 

m,m  m-i,m-i  m-i,m-iv  m,m  m-i,m-i'  '  1 
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and  similarly  for  the  right-hand  plate.  The  equations  may 
now  be  rendered  explicit.  The  solution  is  then  refined  by 
iteration  using  Eqs.3.9  and  3.13.  It  has  been  found  that 
one  iteration  is  generally  sufficient. 

It  is  clear  that  the  above  analysis  is  quite  tedi¬ 
ous,  since  each  case  must  be  analyzed  separately.  In  order 
to  investigate  the  behavior  in  some  more  detail,  the  prob¬ 
lem  was  simplified  as  follows.  Both  plates  were  taken  to  be 
made  of  the  same  material.  Furthermore,  it  was  assumed  that 
no  characteristics  reached  the  interface  from  elsewhere  in 
the  plates  during  the  entire  isentropic  compression  process, 
such  as,  for  example,  characteristics  reflected  from  the  free 
surface  of  the  driver  plate.  The  Euler ian  characteristic 
net  would  therefore  appear  as  in  Fig. 7.1. 

The  problem  can  then  be  rendered  symmetrical  by 
the  transformation 


u  =  U  -  ~  U0 

(7.5a) 

X’  =  R'  -  j  (Ro  +  Ro)  “  \  U0  t 

(7.5b) 

and  similarly  for  the  left-hand  (driver)  plate.  Here  U0  is 
the  initial  driver-plate  velocity, R0  and  R<J  are  the  initial 
surface  positions,  and  initial  conditions  (  )0  apply  at  time 
t  «  0  when  the  driver  plate  motion  starts.  Isentropic  com¬ 
pression  thus  begins  at  the  fourth  shock  reflection,  i.e. 
at  t  =  t4,  X'  =  X4  given  by 

(7.6a) 

£ 

n 

oi 

(7.6b) 
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The  Eulerian  characteristic  net  now  appears  as  in 

Fig. 7.2. 

The  experimental  Hugoniot  (from  Section  IV)  has 
been  fitted  by  the  empirical  expression 


N '  -  *] 


(7.7) 


Thus  the  Riemann  variable  is  given  by  Eq.3.3  as 

c-i 


L  = 


+  1 


w 


(7.8) 


Confining  attention  to  the  right-hand  plate,  we 
have  the  definition 


dX1 

dt 


u 


Introducing  nondimen sional  parameters 


(7.9) 


Eq.7.9  becomes 


(7.10) 


dx  _  u 
cTF  “  Uq 


(7.11) 


Note  that  in  these  coordinates,  impact  in  the  absence  of  air 
would  occur  at  t  «  0,  *  ■  0. 

A  simple  compression  wave  is  induced  in  the  plate 
for  which  the  Riemann  invariant  is 


u  —  L  =  Uq  —  L0  (7.12) 

By  Eqs.7.8  and  7.12,  Eq.7.11  becomes 
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(7.13) 


d*  = 
W 


1  - 


The  surface  pressure  P  is  determined  by  the  pres¬ 
sure  of  the  entrapped  air,  which  is  a  function  of  the  spacing, 
i.e.  by  Eqs.7.1  and  7.2 


(7.14) 


So  that,  finally,  the  differential  equation  of 
motion  of  the  surface  is 

/  ti.  \ 

<y 

1  i  V  . 


dx  _  2  VT 

ar -  1 


¥ 


tt4 


A 


X. 

X 


+  1 


(7.15) 


This  expression  is  best  integrated  numerically. 
Once  %  (t)  has  been  found,  P(t)  may  be  found  directly  from 
Eq . 7 . 14 . 


Independent  parameters  in  the  solution  are 
(u§/V0A)  and  (P0/A),  i.e.,  the  initial  driver  plate  velocity 
and  initial  air  pressure,  nondimens ionalized  with  respect  to 
material  properties  of  the  plates.  Equation  7.7  has  the  in¬ 
teresting  feature  that  most  face-centered  cubic  metals  may 
be  fitted  by  using  a  value  of  %  near  4.3,  while  most  body- 
centered  cubic  metals  may  be  fitted  by  using  a  value  of  £ 
near  2.7.  The  quantities  x4  and  ir4  are,  however,  functions 
of  Uo*  Care  is  necessary  in  scaling  from  one  material  to 
another,  and  such  scaling  is  not  generally  convenient. 
Geometric  scaling,  expressed  by  Eq.7.10  is,  of  course,  very 
useful* 
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7.2  Real  Gas  Analysis 

The  chief  objection  to  Eq.7.15  is  that  air  tempera¬ 
tures  become  very  high,  and  the  perfect  gas  assumption  made 
in  its  derivation  may  be  expected  to  be  poor.  Computations 
based  on  real  gas  behavior  were  therefore  carried  out*  The 
initial  shock  reflection  process  between  rigid  walls  was 

oe 

computed  using  Feldman's  Molier  chart  according  to  the 
method  of  Stollery  and  Maul.^  Results  for  several  cases 
are  given  in  Table  VI  for  comparison  with  results  of  the 
perfect  gas  analysis. 

In  order  to  include  real  gas  behavior  in  the  sub¬ 
sequent  motion,  real  gas  properties  for  equilibrium  air  tabu- 

37  38 

lated  by  Hilsenrath  and  Beckett  *  were  used  in  the  numeri¬ 
cal  integration  in  place  of  Eq.7.14.  In  this  case,  individual 
solutions  must  be  obtained  for  each  initial  pressure,  although 
geometrical  scaling  is  still  valid. 

Solutions  using  the  perfect  gas  analysis  and  the 
real  gas  analysis  are  compared  in  Fig. 7.3,  using  as  an  ex¬ 
ample  aluminum  plates  with  £  *=  4.266,  A  =  188.96  kilobars, 
and  V0  =  0.359  cn^/gm,  and  with  initial  conditions  P0  ■  1.013 
bars,  or  1  atmosphere,  and  U0  =  1.5  mm/nsec.  Pressures  rap¬ 
idly  become  very  high.  For  high  temperature  air  at  very  high 
density,  the  polytropic  gas  law  with  -y  near  1.4  might  be  as¬ 
sumed  to  yield  a  reasonable  approximation  to  the  real  gas 
properties.  A  further  calculation  was  performed,  treating 
the  shock  process  by  using  real  gas  properties,  but  assuming 
the  isentropic  compression  law  to  be  given  by  the  perfect 
gas  law  Eq.7.14.  Thus,  Eq.7.15  was  reintegrated  using  ir4 
and  x.  4  given  by  the  real  gas  analysis  of  Stollery  and  Maul. 
The  resultant  curve  is  also  entered  in  Fig. 7. 3.  It  may  be 
seen  that  while  the  results  of  the  perfect  gas  and  real  gas 
analyses  show  poor  agreement,  the  results  of  the  real  gas 
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shock -perfect  gas  compression  analysis  compare  quite  rea¬ 
sonably  with  those  of  the  real  gas  analysis.  Subsequent 
solutions  were  therefore  obtained  using  real  gas  properties 
for  the  shock  process,  and  perfect  gas  properties  for  the 
subsequent  isentropic  compression. 

Representative  results  are  shown  in  Figs. 7. 4  and 
7.5  in  which  the  initial  driver  plate  velocity  is  varied  at 
a  constant  initial  air  pressure  of  one  atmosphere,  again 
taking  as  an  example  aluminum  plates.  Figure  7.4  shows  sur¬ 
face  trajectories.  In  the  absence  of  air,  the  nondimension- 
alized  trajectory  would  be  a  straight  line  inclined  at  45 
degrees,  the  impact  occurring  at  t  =  0,  =  0.  Figure  7.5 

shows  corresponding  pressure  profiles.  While  the  pressure 
rise  occurs  in  about  the  same  nondimens ional  time  t,  at  each 
of  the  initial  driver  plate  velocities  shown,  the  real  time 
t  is  inversely  proportional  to  the  initial  velocity.  The  ef¬ 
fect  of  the  air  in  degrading  the  wave  profile  is  thus  more 
serious  at  lower  driver  plate  velocities. 

Also  shown  in  Figs. 7. 4  and  7.5  are  profiles  for  an 
initial  driver  plate  velocity  of  1.5  rrmy'n.sec,  but  for  an 
initial  pressure  of  1/100  atmosphere.  As  expected,  the  ef¬ 
fect  of  the  air  is  much  reduced.  It  was  at  one  time  thought 
that  the  effect  of  air  in  an  explosively  driven  plate  experi¬ 
ment,  in  which  the  driver  plate  travel  is  only  about  2  cm, 
could  be  modeled  by  a  plate  experiment  in  an  air  gun,  in 
which  the  driver  plate  travel  is  about  200  cm,  by  reducing 
the  initial  entrapped  air  pressure  in  the  air  gun  in  inverse 
proportion  to  the  initial  spacing.  While  it  is  evident  from 
Eq.7.15  that  such  scaling  is  not  valid,  a  more  graphic  illus¬ 
tration  is  afforded  by  Fig. 7.6  in  which  the  pressure  profiles 
shown  in  Fig. 7. 5  have  been  expanded  to  a  real  time  base. 

(Zero  time  has  been  assigned  arbitrarily  to  each  profile  to 
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allow  presentation  on  a  common  scale.)  The  pressure  profile 
is  considerably  more  degraded  for  the  air  gun  case  (P0  =  1/100 
atm,  2X0  =  200  cm).  Results  of  analyses,  using  real  gas  prop¬ 
erties  throughout,  are  also  shown  for  comparison. 

7.3  Subsequent  Wave  Profiles 

Thus  far,  only  pressure  profiles  at  the  interface 
have  been  considered.  A  gradually  increasing  compression 
wave  is  induced  in  both  the  target  and  driver  plates.  As 
the  compression  waves  propagate  into  the  interior  of  the 
plates,  they  will  steepen  and  eventually  form  shock  waves. 
Thus,  at  some  distance  in  the  interior,  the  effect  of  the 
air  is  no  longer  discernible.  In  order  to  follow  this  pro¬ 
cess,  the  waves  propagating  into  the  plates  were  analyzed  by 
the  method  of  characteristics  outlined  in  Section  IV,  treating 
shock  waves  as  isentropic  compressions  in  the  weak  shock  ap¬ 
proximation.  Resultant  pressure-real  time  profiles  at  vari¬ 
ous  stations  in  the  plates  are  shown  in  Fig. 7. 7  for  each  of 
the  four  initial  velocities  investigated.  Thus,  at  an  initi¬ 
al  driver  plate  velocity  of  1.5  rar0isec  for  example,  the 
compression  wave  has  steepened  into  a  shock  in  2  mm  of 
travel,  except  for  a  small  forerunner.  However,  at  an 
initial  driver  velocity  of  0.35  ran/nsec,  a  shock  is  not 
fully  formed  until  the  pulse  has  propagated  40  mm. 

If  the  driver  plate  thickness  is  greater  than  the 
distance  which  the  compression  wave  must  travel  in  order  to 
steepen  into  a  shock,  then  no  effect  due  to  the  air  will  be 
apparent  on  the  subsequent  wave  profile  travelling  into  the 
target  (except  for  a  small  effect  due  to  the  forerunner) . 
However,  if  the  driver  plate  thickness  is  less  than  this 
distance,  then  the  release  wave  reflected  from  the  free  sur¬ 
face  of  the  driver  plate  will  be  affected.  Thus,  if  the 
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compression  wave  arriving  at  the  free  surface  of  the  driver 
plate  has  some  profile  of  finite  width,  this  profile  is 
transferred  to  the  release  wave  reflected  back  into  the 
driver  plate.  The  release  wave  now  will  no  longer  steepen, 
but  degrade  further,  and  on  passing  into  the  target  will  be 
expected  to  show  a  different  profile  from  that  induced  in 
the  target  in  the  absence  of  air. 

An  analytical  solution  of  the  problem  requires  so¬ 
lution  of  the  characteristic  net  resulting  from  the  reflec¬ 
tion  of  the  compression  wave  at  the  driver  plate  free  sur¬ 
face.  A  sample  calculation  has  been  carried  through  for  the 
following  initial  conditions  P0  =  1  atm,  2X0  -  2  cm, 

U0  =  1.5  mn/psec,  and  driver  plate  thickness  0.8  mm, again 
with  aluminum  plates.  The  characteristic  mesh  is  shown  in 
Fig. 7. 8. 

The  corresponding  pressure-real  time  profile  at 
the  interface  is  shown  in  Fig. 7.9.  The  solution  (assuming 
an  isentropic  shock)  for  the  no-air  case  is  also  shown  for 
comparison.  The  pressure-Eulerian  distance  profile  is  shown 
in  Fig. 7. 10,  at  a  time  0.357  microseconds  after  impact  in 
the  absence  of  air  (13.693  microseconds  after  start  of 
driver  plate  motion) .  It  may  be  seen  that  the  only  signifi¬ 
cant  difference  introduced  by  the  presence  of  air  is  a 
rounding  of  the  head  of  the  rarefaction  wave.  However,  this 
rounding  leads  to  important  differences  as  the  rarefaction 
catches  up  to  the  shock,  and  the  shock  is  attenuated.  In 
particular,  the  peak  pressure  will  be  attenuated  much  sooner, 
and  the  peak  pressure-distance  curve  will  not  exhibit  a  well- 
defined  "catch-up"  distance,  the  comer  (see  Section  V)  be¬ 
ing  rounded. 

The  passage  of  the  release  wave  across  the  inter¬ 
face  containing  the  air  layer  was  accomplished  in  the  above 
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calculation  by  simply  ignoring  the  presence  of  the  air,  and 
allowing  the  characteristics  to  propagate  across  the  inter¬ 
face  without  refraction.  While  this  is  a  good  approximation 
during  the  initial  part  of  the  pressure  release  process,  it 
is  not  strictly  valid  for  the  final  part.  As  the  pressure 
at  the  interface  drops,  the  air  will  expand,  causing  the 
driver  plate  to  rebound.  A  partial  analysis  was  undertaken 
in  which  the  correct  interface  conditions  (Eqs.3.8,  3.9, 

3.10,  3.11,  3.13,  7.1,  and  7.2)  were  used,  iterating  at  each 
step.  The  solution  followed  the  profile  of  Fig, 7. 9  down  to 
about  10  kb,  at  which  point  the  divergence  became  serious. 
The  characteristic  mesh  is  quite  complex  due  to  further  wave 
reflections  from  the  driver  free  surface,  and  the  solution 
was  discontinued  at  this  point. 

In  the  previous  calculation,  the  first  release 
wave  from  the  free  surface  of  the  driver  plate  arrived  back 
at  the  interface  after  the  compression  of  the  air  had  been 
completed,  (see  Fig. 7.8).  Thus,  the  compression  was  handled 
by  Eq,7.15  and  the  reflection  of  the  compression  wave  at  the 
driver  free  surface  was  handled  separately.  More  serious 
degradation  of  the  wave  profile  is  to  be  expected  if  the  re¬ 
flected  release  wave  arrives  at  the  interface  before  the 
peak  pressure  has  been  attained  there.  Equation  7.15  can 
now  no  longer  be  used,  and  a  characteristic  method  using  the 
correct  conditions  at  the  interface  is  required.  A  wave  pro¬ 
file,  computed  on  this  basis  for  the  following  initial  con¬ 
ditions  is  shown  in  Fig. 7.11;  P0  =  1  atm,  U0  ■  0.35  nm^M-sec, 
2X0  ■*  2  cm,  driver  plate  thickness  1  mm.  The  peak  pressure 
attained  is  22  per  cent  below  that  attained  in  the  absence 
of  air,  or  with  a  much  thicker  driver  plate. 

There  will  thus  be  three  types  of  behavior  which 
might  be  observed.  If  the  driver  plate  is  sufficiently 
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thick,  virtually  no  effect  due  to  the  air  will  be  apparent, 
since  the  compression  steepens  into  a  shock  before  reaching 
the  free  surface  of  the  driver  plate.  For  thinner  driver 
plates,  the  wave  induced  in  the  target  will  show  some 
rounding  at  the  head  of  the  rarefaction,  and  hence  will  at¬ 
tenuate  somewhat  more  rapidly  as  the  rarefaction  overtakes 
the  shock  wave.  For  very  thin  targets,  strong  rarefaction 
waves  arrive  at  the  interface  from  the  rear  surface  of  the 
driver  plate  before  the  peak  pressure  has  been  attained, 
and  the  pressure  attained  in  the  absence  of  air  is  never  at¬ 
tained  at  the  interface.  Approximate  boundaries  can  be 
found  for  the  regions  in  which  each  of  these  types  of  be¬ 
havior  may  be  expected,  based  on  the  previous  calculations. 
These  boundaries  are  shown  in  Fig. 7. 12  for  the  range  of  ve¬ 
locities  investigated  in  the  present  examples,  and  for 
aluminum  driver  plates,  with  an  initial  air  pressure  of 
1  atmosphere.  Geometric  scaling  is  valid,  and  the  driver 
plate  thickness  has  thus  been  nondimensionalized  with  respect 
to  one  half  the  initial  plate  spacing  (X0).  Points  repre¬ 
senting  the  conditions  corresponding  to  the  profiles  of  Figs. 
7.9  and  7.11  have  been  inserted. 

It  must  be  noted  that  this  diagram  cannot  be  scaled 
to  other  materials  or  initial  pressures,  and  in  general,  it 
is  highly  advisable  to  carry  out  a  complete  calculation  for 
each  case  of  interest. 

In  order  to  assess  whether  the  effect  of  side 
leakage  is  important  near  the  center  of  the  plates,  calcula¬ 
tions  were  made  of  the  distance  travelled  by  sound  waves 
during  the  compression  cycle.  Thus  the  quantity 
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where  c  “  U  'i's  son*c  velocity  in  the  air,  was  summed 
over  all  time  steps  n  in  the  numerical  integration  of  Eq. 
7.15,  from  time  t4 .  For  each  of  the  four  initial  velocities 
considered  above  at  an  initial  pressure  of  one  atmosphere, 
the  nondimens ional  distance  5  was  in  excess  of  10.  Thus, 
for  example,  for  plates  5  x  5  cm  in  lateral  extent,  with  a 
2  cm  initial  spacing,  release  waves  from  the  edges  are  able 
to  travel  to  the  center  and  back  at  least  twice  during  the 
compression  cycle.  The  idealized  one- dimens ional  analysis 
therefore  provides  only  an  upper  limit.  In  fact,  consider¬ 
able  leakage  must  occur  in  most  practical  cases.  However, 
without  a  detailed  three-dimensional  analysis  of  the  air 
flow,  it  is  not  possible  to  determine  the  quantitative  de¬ 
parture  from  the  results  of  the  idealized  one- dimens ional 
analysis. 

7.4  Summary 

It  is  convenient  to  summarize  the  main  conclusions 
of  the  analysis  as  follows: 

1.  In  analyzing  the  air  effect,  real  gas  effects 
must  be  included.  Assuming  perfect  gas  behavior  for  the 
shock  compression  phase  leads  to  large  errors. 

2.  The  phenomenon  may  be  scaled  geometrically, 
in  which  case  time  must  be  scaled  with  the  geometric  scale 
ratio.  Scaling  of  initial  pressure  and  velocity  involves 
plate  material  properties.  Initial  pressure  and  velocity, 
nondimens ional ized  by  material  properties,(P0/A)  and 
(uo/V0A) ,are  independent  parameters. 

3.  Degradation  of  wave  profiles  is  more  serious 
for  lower  initial  velocities,  and  for  more  compressible 
materials  (lower  A).  While  the  effects  should  be  quite 
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small  for  most  normal  metal  plate  systems,  except  for  very 
thin  driver  plates,  the  effect  may  be  much  more  important 
for  plastics. 

4,  Magnitude  of  the  effect  depends  on  the  driver 
plate  thickness,  other  initial  parameters  being  constant. 

For  relatively  thick  driver  plates,  no  effect  will  be  ob¬ 
servable  once  the  entire  pulse  has  entered  the  target.  For 
thinner  driver  plates,  the  only  important  effect  will  be 
the  rounding  of  the  head  of  the  release  wave  entering  the 
target,  with  a  consequent  large  reduction  in  catch-up  dis¬ 
tance  (start  of  attenuation  of  the  shock  by  the  release 
wave)  and  subsequent  more  rapid  attenuation  of  the  shock  by 
the  release  wave.  For  very  thin  driver  plates,  the  peak 
pressure  at  the  interface  will  be  attenuated  by  release 
waves  reflected  from  the  rear  surface  of  the  driver  plate, 
and  the  peak  pressure  induced  in  the  absence  of  air  will 
never  be  attained.  As  the  pulse  propagates  into  the  target, 
subsequent  attenuation  is  also  very  much  more  rapid. 

5,  Since  the  most  important  effect  of  air  is  in 
altering  the  wave  profile,  it  is  not  possible  to  counteract 
the  effect  by  increasing  the  driver  plate  velocity.  Such 
action  will  only  lead  to  a  distorted  wave  of  higher  ampli¬ 
tude,  which  cannot  have  the  same  attenuation  characteristics 
as  a  pulse  induced  in  the  absence  of  air. 

6,  In  almost  all  practical  cases,  it  is  expected 
that  side  leakage  of  the  air  will  be  important,  and  there¬ 
fore  the  foregoing  one- dimensional  analysis  will  provide 
only  an  upper  limit  on  the  effect  near  the  center  of  the 
plates.  Leakage  will  affect  areas  near  the  edges  of  the 
plates  to  a  much  greater  degree  than  in  the  center,  and  the 
wave  profile  will  more  nearly  approach  that  induced  in  the 
absence  of  air.  The  result  will  be  a  complex  three-dimen- 
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sional  nonplanar  wave  system  propagating  into  the  target. 
Nonplanarity  of  the  wave  system  may  in  itself  affect  such 
phenomena  as  spall. 

7.  In  the  light  of  the  foregoing,  it  is  concluded 
that  if  the  present  analysis  shows  that  the  air  effect  is 
important  in  a  particular  series  of  tests,  then  the  resultant 
data  from  those  tests  is  unreliable.  The  only  course  is  to 
remove  the  air  in  the  experiment. 
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SECTION  VIII 


CONCLUSIONS 

In  an  analysis  of  the  effect  of  a  distributed  in¬ 
tense  impulsive  load  on  a  structure,  it  is  necessary  to 
compute  the  stress  wave  propagation  through  the  thickness  of 
the  structure.  Two  methods,  which  have  been  discussed,  are 
available  viz.  the  finite  difference  method  (Section  II )  and 
the  characteristics  method  (Section  III).  Furthermore 
several  approximate  constitutive  equations  have  been  con¬ 
structed  (Section  IV)  which  are  valid  in  different  ranges. 

For  two  aluminum  alloys  for  which  data  exists, 
little  or  no  increase  in  yield  strength  is  apparent,  and  it 
has  been  possible  to  relate  the  dynamic  constitutive  equa¬ 
tion  for  moderate  stresses  to  static  material  properties  via 
second  order  finite  strain  theory. 

Success  of  the  analysis  of  stress  wave  propagation 
depends  on  the  correct  choice  of  computation  method  and  con¬ 
stitutive  equation.  Several  examples  have  been  given 
(Section  V) ,  but  no  general  rules  can  be  formulated. 

Spall  fracture  has  been  discussed  in  some  detail 
(Section  VI) .  For  copper,  for  which  some  unambiguous  experi¬ 
mental  data  is  available,  only  a  moderate  increase  in  frac¬ 
ture  threshold  is  observed  at  the  high  rates  involved  in 
spallation.  The  fracture  threshold  must  be  determined  ex¬ 
perimentally,  the  plate  impact  experiment  being  ideally 
suited  to  this  purpose. 

The  effect  of  air  in  the  plate  impact  experiment 
has  been  analyzed  in  detail  (Section  VII),  but  it  is  only 
possible  to  obtain  an  upper  limit  on  the  effect.  Therefore, 
if  an  analysis  shows  the  air  eff°^t  to  be  important  in  a 
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particular  experiment,  the  analysis  cannot  be  used  to  cor¬ 
rect  the  resultant  data,  rather  the  experiment  should  be  re¬ 
peated  in  a  vacuum. 

For  materials  other  than  the  aluminum  and  copper 
alloys  for  which  data  exists,  further  experimental  work  is 
required  to  determine  if  their  behavior  is  similar.  The 
constants  in  the  constitutive  equation  and  fracture  criteria 
for  these  materials  must  be  evaluated  experimentally  before 
an  analysis  of  stress  wave  propagation  and  spallation  can  be 
undertaken. 
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►  r 

Fig. 3.1  Characteristic  Net  for  Reflection  at  an  Interface 


Fig. 3. 2  Characteristic  Net  for  Formation  and  Reflection  of 
a  Weak  Shock 
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Fig. 4 .2a  Offset  of  the  Hugoniot  from  the  Adiabat  In  Pure 

Dilatation  for  2024  Aluminum  (data  from  Fowles^) 
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Fig. A. 2b  Offset  of  the  Hugoniot  from  the  Adiabat  in  Pure  Dilatation 

18 

for  6061-T6  Aluminum  (data  from  Lundergan  ) 
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Fig. A. 3  Murnaghan  Fit  to  the  Dynamic  Adiabat  for  Aluminum 
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Fig. 5.1  Lagrangian  Characteristic  Mesh  for  Plate  Impact 
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ig.5.2a.  Wave  Profiles  at  2  [isec  After  Impact.  Copper 
Driver  and  Target  Plates  with  an  Impact  Ve¬ 
locity  of  12.6  x  10s  cn/sec 
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Fig. 5. 2b.  Wave  Profiles  Continued.  At  3  nsec  After  Impact 
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Fig. 5. 2c.  Wave  Profiles  Continued.  At  4  nsec  After  Impact 
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Fig. 5. 2d.  Wave  Profiles  Continued.  At  5  M-sec  After  Impact 
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Wave  Profiles  Continued,  At  6  M-sec  After  Impact 
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Fig.5.2f.  Wave  Profiles  Concluded.  At  7  ixsec  After  Impact 
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Fig. 5. 3a.  Comparison  of  Characteristic  and  Finite  Dif¬ 
ference  Methods  for  Plate  Impact  Using  the 
Murnaghan  Equation 
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Fig. 5. 3b.  Concluded 
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Fig. 6.1  Sections  Through  1/4-inch-thick  E.T.P  Half -Hard  Copper 
Target  Plates  Struck  by  1/8- inch- thick  Driver  Plates 
of  the  Same  Material.  Ferric  Chloride  Etch 


Fig. 6. 2  Tension  Test  Arrangement 
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Driver  Plate  Velocity  (10*  cny'sec) 

Fig. 6. 4  Spall  Velocity  as  a  Function  of  Driver  Plate  Velocity 
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Fig. 6. 5  Characteristic  Mesh  for  the  Elastic-Plastic 
Wave  System  in  the  Sandia  Plate  Experiments 
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Fig. 6. 6a  Stress-Time  History  at  a  Station  at  the  Target  Centre 
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Fig. 6. 6b  Stress-Time  History  at  a  Station  One  Quarter  of  the 
Target  Thickness  from  the  Interface 
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Fracture  Time  tfi 

Fig. 6. 7  Qualitative  Upper  Fracture  Threshold  Time  Dependence  for 
Half -Hard  E.T.P. Copper 


Fig. 6. 8  Edge  Effect  With  Groove 
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Fig.  7.1  Characteristic  Mesh  for  Plate  Impact 
with  Air  Present 
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Fig.  7.2  Non-dimens ionalised  Characteristic  Mesh 
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Fig. 7.3  Comparison  of  Perfect  Gas  and  Real  Gas  Analyses  of  the  Air  Effect 


Fig. 7.4  Surface  Trajectories  in  the  Presence  of  Air 
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Fig. 7. 5  Interface  Pressure  Profiles  in  the  Presence  of  Air 
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Fig.  7.6  Effect  of  Air  Between  Impacting  Plates 
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Fig. 7. 7a  Wave  Profiles  in  the  Target 
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Fig. 7. 7c  Wave  Profiles  in  the  Target 
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Fig. 7.8  Eulerian  Characteristic  Mesh  for  Impact  with  Air 
Between  Driver  and  Target 
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Fig.  7.9  Wave  Profile  at  Interface  Relative  to  No -Air  Impact  Time 
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Fig.  7.10  Wave  Profiles  at  Fixed  Time  After  No-Alr  Impact 
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Fig. 7. 11  Pressure  Rise  at  the  Interface  for 
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Fig.  7.12  Approximate  Regions  of  Effects  of  Air  in  Piste  Impact 
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Dnlitdrlt  Principal  Cmmwia 

i 

: 

i 

* 

i1 

i 

■ 

! 

; 

11 attic 

dQa  1  fr?"B  doa  *  S  *  S  da* 

affi-g  da* 

11*1  m  do«  -  dfl»  *  aVfe  *>a 

Plaatlc* 

-  v  *  *  * 

|  ?  it  ;  Jffj  il 

^  ^  dt 

+  §  dT  i  ¥  j  d*  :  T  ^  JT 

1 

a 

£ 

2 

tl- 

14 

lluElC 

dij  !  tj  :  o 

4-“.  !  -4*.  >  -4*. 

Platt 1c j 

Klaatlc  C«pPn«at 

*S  =  i*I  *  i*i 

^  (Ida;  -  <ta*i  j  -  |  <*<*;  -  dej]  |  -  1  <2<n£  - 

riaatlc ; 

r  late  lc  Caaponanc 

«:  f  -4*:  <  -4«: 

0 

*;  i  -4*;  i  -4*; 

*Tbt  L±pp*r  alp  la  tppfoprlu*  i**n  Aoa  la  poiltlm  tba  lowar  *L*n  w*vp»  dki^  la  iHtatLva, 
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TABLE  VI 


Multiple  Shock  Reflection  Parameters 

(C.G.S.  units  throughout) 


PERFECT  GAS 

REAL  GAS 

Vo 

Region 

P 

P 

T 

P 

P 

T 

xiO4 

xlO6 

xlO*3 

°K 

xlO6 

xlO"3 

°K 

15 

0 

1.013 

1.173 

300 

* 

* 

* 

1 

34.54 

5.998 

2000 

* 

* 

* 

2 

233.4 

19.53 

4151 

213.8 

25.70 

3300 

3 

971.0 

49.91 

6760 

902.9 

70.78 

4650 

4 

3084.9 

109.1 

9820 

2796.0 

195.0 

5900 

10 

0 

1.013 

1.173 

300 

* 

* 

* 

1 

16.17 

5.169 

1086 

* 

* 

* 

2 

93.27 

15.640 

2071 

* 

* 

* 

3 

357.8 

38.19 

3255 

327.8 

40, 26 

2600 

4 

1080.1 

80.93 

4636 

978.6 

83.16 

3800 

7 

0 

1.013 

1.173 

300 

* 

* 

* 

1 

8.931 

4.267 

727 

* 

* 

* 

2 

41.91 

11.63 

1251 

* 

* 

* 

3 

143.2 

26.57 

1872 

* 

* 

* 

4 

400.6 

53.70 

2591 

359.4 

57.53 

2250 

3-5 

0 

1.013 

1.173 

300 

* 

* 

* 

1 

3.541 

2.714 

453 

* 

* 

* 

2 

10.05 

5.538 

630 

* 

* 

* 

3 

24.70 

10.31 

832 

★ 

* 

* 

4 

54.48 

17.88 

IO58 

★ 

* 

* 

*  Same  as  perfect  gas  conditions. 
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